# Geometry, percolation and transport properties of fracture networks derived from line data

S. Sisavath, V. Mourzenko, P. Genthon, J. -F. Thovert, P. M. Adler

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

► To cite this version:

S. Sisavath, V. Mourzenko, P. Genthon, J. -F. Thovert, P. M. Adler. Geometry, percolation and transport properties of fracture networks derived from line data. Geophysical Journal International, 2004, 157, pp.917-934. 10.1111/j.1365-246X.2004.02185.x. insu-03600286

HAL Id: insu-03600286

https://insu.hal.science/insu-03600286

Submitted on 7 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_2.png" /></figure>
# Geometry, percolation and transport properties of fracture networks derived from line data

S. Sisavath,<fnref n="1" /> V. Mourzenko,<fnref n="2" /> P. Genthon,<fnref n="3" /> J.-F. Thovert<fnref n="2" /> and P. M. Adler<fnref n="1" />

Accepted 2003 October 21. Received 2003 October 6; in original form 2002 October 14

## SUMMARY

In most geological instances, 2-D or 3-D fracture distributions are not available from field data. We show here that when data relative to fractures are collected along a line such as a road or a well, estimations can be given to the major geometrical properties of the corresponding fracture networks, such as the volumetric density of fractures, their percolation character and their macroscopic permeability. All these formulae are analytical and can be split into two parts; the first one can be derived from the measured data, while the second one requires some assumption on the lateral extension of the fractures and on their permeability. All these techniques are applied to fractures located in the Baget watershed. They are also validated on a granite block whose structure is fully known. Extensions are proposed for networks with variable permeabilities and polydisperse fractures.

**Key words:** fracture network, geometry, line data analysis, percolation, permeability, reconstruction.

## 1 INTRODUCTION

The macroscopic properties of a fracture network consist of its geometrical and topological properties (such as connectivity and percolating character), and of its transport properties (such as permeability).

The main objective of this work is to show how a fracture network can be characterized by data collected along a line. This line can be located at the surface of the ground, but it can also be a well. Networks can be subsequently reconstructed according to these data and the geometrical properties, namely connectivity and percolation, can be estimated; moreover, transport properties such as permeability can be determined as well. In other words, we wish to show that relatively important overall information can be extracted from relatively poor data.

Fractures are usually generated in a random way (e.g. Koudina et al. 1998), whereas their positions and characteristics are likely to depend on the geological features of the region under study. Another technique is stereological analysis which was used for instance by Berkowitz & Adler (1998); the input data were traces collected on a plane surface. In this work, we focus on fractures located in the Baget watershed (in the southeast of France) which presents the peculiarity of being karstified so that an analysis of the fracture distribution could provide some hints of the drainage pattern of the basin (Pistre et al. 1999); this means that the transport properties of the fracture networks will be dependent, to a great extent, on phenomena such as deposition and dissolution. Due to bad outcropping conditions, fractures are only visible along recently excavated

roads; hence, geometrical parameters of the fractures such as their dips and strikes were mapped along these roads. These parameters were then included into our network model. It should be stressed that the karstified character will be studied in a different contribution.

This paper is organized as follows. Section 2 presents the major theoretical considerations. After a few definitions, the volumetric density of an event is determined as a function of measurable quantities. Then, it is shown how the networks can be reconstructed when the individual events are assumed to be randomly distributed in space. The general percolation properties of these networks can be derived from quantities such as the mean number of intersections per event, which can be partly expressed in terms of the measured data. The determination of the permeability of these networks is presented and two general estimations are proposed: the first one derived from percolation theory is valid close to the percolation threshold, while the second one is valid for high densities. It should be emphasized that many of these properties can be derived from relatively straightforward analytical formulae whose simplicity is a decisive factor for the applications.

Section 3 applies these concepts to the particular case of the Baget watershed. After a detailed analysis of the line surveys, several types of fracture networks are reconstructed; fractures may be of the same size (monodisperse fractures) or of two different sizes (bidisperse fractures), and each event can be composed of one or several parallel fractures. Examples of networks are presented. Their percolation properties are discussed as well as their flow properties.

The first part of Section 4 is devoted to the validation of the methodology on a block of granite whose structure and properties

<footnote n="1">¹IPGP, 4 Place Jussieu, 75252-Paris Cedex 05, France. E-mail: adler@ipgp.jussieu.fr</footnote>

<footnote n="2">²LCD, SP2MI, BP 179, 86960 Futuroscope Cedex, France</footnote>

<footnote n="3">³Observatoire Midi-Pyrénées, 14 Avenue E. Belin, 31400-Toulouse, France</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Notations: (a) intersection of the observation line L with a series of fractures; (b) intersection of a fracture with an observation line and the corresponding intersection volume.</figcaption></figure>

are known extensively (Ledésert et al. 1993; Gonzalez Garcia et al. 2000). In the second part of Section 4, it is first recalled how the permeability of the individual fractures can be determined; the approach and the formulae presented in Section 2 are extended to networks with variable permeabilities, and applications to the Baget watershed are detailed. Finally, the polydisperse character of the fracture sizes is addressed; an extension of the classical Snow formula is proposed.

Some concluding remarks end this paper.

## 2 THEORY

### 2.1 Definitions

Let us consider a line L (cf. Fig. 1a) which is crossed by fractures; the line is parallel to the unit vector **p**. More generally, since these features include single fractures, series of parallel fractures and possibly faulted zones, these intersections are classified as events numbered by *f* for which the position $x_f$ along the profile and the width $w_f$ can be measured; $w_f$ is of course equal to zero when the event contains a single fracture. In the opposite case, the event *f* contains $n_{f,p}$ fractures which are more or less parallel with an average spacing $d_f$ with the obvious relation

$$w_f = n_{f,p} d_f. \qquad (1)$$

In addition to this, the orientation of each fracture should be measured. Instead of the strike and dip angles, we use directly the unit vector $\mathbf{n}_i$ perpendicular to the fracture plane, which can be deduced from the information above; the vertical component of these unit vectors is taken as positive.

Another useful quantity is the spacing $s_i$ between two successive events $i$ and $i + 1$ along L

$$s_i = x_{i+1} - x_i. \qquad (2)$$

### 2.2 Spatial distributions and volumetric density

In the subsequent simulations, we aim at generating random realizations of fracture networks in agreement with the field observations. To this end, it is necessary to know how many events of each type are to be inserted in a given sample volume.

To determine this from the line surveys, we need to know the density of a given event characterized by the subscript *f*. Suppose that this event is a fracture of area $A_f$ and of normal $\mathbf{n}_f$. This fracture intersects a segment of length *L* on a line parallel to **p** if its center

belongs to the cylinder of base $A_f$ and length *L* as illustrated in Fig. 1(b). The volume of this cylinder is

$$V = L A_f |\mathbf{p} \cdot \mathbf{n}_f|. \qquad (3)$$

If $\rho_f$ is the number of such events per unit volume, the number of intersections $n_f$ of this type of event with *L* is given by

$$n_f = \rho_f L A_f |\mathbf{p} \cdot \mathbf{n}_f|. \qquad (4)$$

Therefore, the average spacing $s_f$ between two intersections is given by

$$s_f = \frac{L}{n_f} = \frac{1}{\rho_f A_f |\mathbf{p} \cdot \mathbf{n}_f|}. \qquad (5)$$

Actually, this quantity $s_f$ is measurable along the profile *L*. It can be used to express the volumetric density of the event *f* as

$$\rho_f = \frac{1}{L_f A_f} \text{with } L_f = s_f |\mathbf{p} \cdot \mathbf{n}_f|. \qquad (6)$$

This apparently formal decomposition has the great advantage of dividing $\rho_f$ into two terms, namely $L_f$ which is known from the measurements along the line and $A_f$ whose value has to be hypothesized. Such a decomposition will be used repeatedly in this paper, and it can be considered as one of the major contributions of this work.

Note that according to (5), $1/L_f$ is also equal to $\rho_f A_f$, i.e. to the area of the event *f* per unit volume.

The main problem for the determination of $\rho_f$ is the estimation of the area $A_f$. Here, it will remain an unknown quantity since data collected along a line cannot give any information on their extent. Note that when 2-D trace maps are available, these quantities can be estimated after a stereological analysis (Berkowitz & Adler 1998).

### 2.3 Numerical reconstruction

It is useful to give some general information on the reconstruction procedure which involves two steps.

First, the events are assumed to be Poissonian (a hypothesis subject to verification) and they are inserted at random locations: their number results from the volumetric density $\rho_f$ in (6) and from the sample volume. Secondly, whenever it is relevant, each event is decomposed into series of parallel fractures, the number and spacing of which are deduced from the data *w* and *d*.

The first step of the generation follows the procedure of Huseby et al. (1997). The generated samples are made up of cubic unit cells of size *L*. In order to minimize size effects in the determination of the network transport properties, periodicity conditions are applied,
so that an infinite medium results from the juxtaposition of identical unit cells.

The events and the fractures they contain are assumed to have identical plane circular or polygonal shapes. Their size is quantified by the radius R of the circumscribed circle.

## 2.4 Connectivity

Within the framework of Poissonian spatial distribution of the fracture events, much can be deduced about the connectivity of the network from the line survey data.

Recall firstly the definition of the exclusion volume $V_{ex}(1, 2)$ of two objects $f_1$ and $f_2$ (Balberg et al. 1984). It is the volume surrounding $f_1$ in which the center of $f_2$ must be located in order for them to intersect. Note that the definition is symmetric, with $V_{ex}(1, 2) = V_{ex}(2, 1)$. Thus, if $\rho_i$ is the volumetric density of objects $f_i$, $\rho_i V_{ex}(i, j)$ is the mean number of $f_i$ objects intersecting $f_j$. By simple summation, the total number $\rho_I$ of object intersections per unit volume is

$$ \rho_I = \frac{1}{2} \sum_i \sum_j \rho_i \rho_j V_{ex}(i, j). \qquad (7) $$

For the sake of simplicity, we suppose that the fracture events can be regarded as plane convex polygons, i.e. that their width $w_f$ is small compared to their extension. The argument can be easily generalized to events with non-zero thickness, by modifying the expression (8) below for the excluded volume.

It can be shown (Adler & Thovert 1999) that if the objects $f_i$ are plane convex polygons, with areas $A_i$, perimeters $P_i$, and relative orientations given by the angle $\gamma_{i,j}$ between their normal vectors, the excluded volume is given by

$$ V_{ex}(i, j) = \frac{\sin \gamma_{i,j}}{\pi} (A_i P_j + A_j P_i) \text{ with } \cos \gamma_{i,j} = \mathbf{n}_i \cdot \mathbf{n}_j. \qquad (8) $$

By injecting (8) into (7), the density of intersections reads

$$ \rho_I = \frac{1}{2} \sum_i \sum_j \rho_i \rho_j \frac{\sin \gamma_{i,j}}{\pi} (A_i P_j + A_j P_i). \qquad (9) $$

However, $\rho_i$ is directly related to the volumetric area $1/L_i$ by eq. (6). Hence,

$$ \rho_I = \frac{1}{2\pi} \sum_i \sum_j \frac{\sin \gamma_{i,j}}{L_i L_j} \left( \frac{P_i}{A_i} + \frac{P_j}{A_j} \right). \qquad (10) $$

If the shapes and sizes of the various events are known, they can be used in the fully general expression (10). Partial summations can also provide more detailed statistics, such as the mean numbers of intersections between specific events or families of events.

In the absence of such information, we suppose here that all the fracture events have the same shape and size, and thus, the same area A and perimeter P. The density of intersections then reads

$$ \rho_I = \frac{P}{\pi A} \sum_i \sum_j \frac{\sin \gamma_{i,j}}{L_i L_j}. \qquad (11) $$

On the other hand, the total volumetric density $\rho$ of events, i.e. the number of fractures per unit volume, is simply the sum of the densities $\rho_i$:

$$ \rho = \frac{1}{A} \sum_i \frac{1}{L_i}. \qquad (12) $$

Moreover, the mean number $\rho'$ of intersections per event, which directly quantifies the connectivity of the network and will be used hereafter as a dimensionless measure of the network density is

$$ \rho' = \frac{2\rho_I}{\rho} = \frac{2P}{\pi} \frac{\sum_i \sum_j \frac{\sin \gamma_{i,j}}{L_i L_j}}{\sum_i \frac{1}{L_i}}, \qquad (13) $$

where the factor 2 comes from the fact that an intersection necessarily belongs to two different fractures.

A global effective exclusion volume $V_{ex}$ can be defined from $\rho$ and $\rho'$, as

$$ V_{ex} = \frac{\rho'}{\rho} = \frac{2}{\pi} AP \frac{\sum_i \sum_j \frac{\sin \gamma_{i,j}}{L_i L_j}}{\left(\sum_i \frac{1}{L_i}\right)^2}. \qquad (14) $$

For an isotropic distribution of monodisperse objects, $V_{ex} = AP/2$ (Adler & Thovert 1999).

As (6), all these formulae can be decomposed into a measurable part and a part which has to be hypothesized.

## 2.5 Probability of percolation

It was shown in Section 2.4 that the connectivity depends both on the density and size of the fractures. However, percolation is a topological property, and it is therefore desirable to quantify the event density by an intrinsic, topological, dimensionless parameter.

The mean number of intersections $\rho'$ defined by (13) is the ideal quantity in this respect. It can be regarded both as a topological parameter, which directly measures the connectivity of the network, and as a dimensionless density, since it is the mean number of objects per volume $V_{ex}$ (see eq. 14). It has been successfully used by Huseby et al. (1997) to describe and unify the topological and geometrical properties of isotropic fracture networks, and by Koudina et al. (1998) for the flow properties of such networks. In the following numerical results, it is directly measured on the reconstructed samples as the ratio of the actual number of fracture intersections $N_I$ to the number of fractures $N_f$

$$ \rho' = \frac{2N_I}{N_f}. \qquad (15) $$

The factor 2 is due to the fact that an intersection belongs to two fractures.

Percolation is defined as the existence of a spanning continuous path, where a fluid can circulate, across the medium. Percolation is a crucial topological property which controls many geometrical and transport properties of the network. This concept originated and was closely studied in discrete sites or bonds lattices (see, e.g. Stauffer & Aharony 1994), where the density of occupied (or conducting) sites or bonds is measured by a probability (or concentration) p. In large systems, percolation occurs when p exceeds a critical value $p_c$, known as the percolation threshold. For p close to $p_c$, many geometrical or transport coefficients are known to scale as power laws of the difference $p - p_c$, according to the standard form

$$ X \propto (p - p_c)^a. \qquad (16) $$

The quantity X may represent the correlation length, the fraction of sites connected to the infinite cluster, or the transport coefficients of the system. Different exponents are associated with the various quantities, but each is generally believed to be universal, i.e. insensitive to the details of the underlying lattice.

It is one of the great merits of the concept of excluded volume and of the dimensionless density $\rho'$, to allow for a transposition of the
on-lattice percolation theory in terms of p to continuous percolation systems, with a formulation in terms of $\rho'$. A survey of this topic is provided by Adler & Thovert (1999). In particular, a critical density $\rho'_c$ exists for a given class of continuous systems, which plays the role of $p_c$ in lattice systems. For networks of randomly-oriented fractures, Huseby et al. (1997) determined this percolation threshold to be

$$ \rho'_c \sim 2.26. \qquad (17) $$

Of course, this threshold is different for non-isotropic networks. However, it provides an order of magnitude to which the estimation (13) can be compared; it can be easily decided a priori if a given fracture network is likely to percolate or not.

## 2.6 Flow properties

The flow properties of these networks can be studied by the methods and numerical tools described by Koudina et al. (1998). The solid matrix containing the fractures is assumed to be impervious.

Let us first consider a single fracture at a local scale characterized by a typical aperture $b_0$, which is assumed to be much smaller than the lateral extent R of the fracture. Then, the flow of a Newtonian fluid at low Reynolds number is governed by the Stokes equations.

At a scale $\mathcal{L}$, which is intermediate between $b_0$ and R, the flow is governed by the Darcy equation

$$ j = -\frac{1}{\mu} \sigma \cdot \overline{\nabla p}, \qquad (18) $$

where j and $\overline{\nabla p}$ are the locally averaged flow rate per unit width [$L^2T^{-1}$], and the pressure gradient $\sigma[L^3]$ is the fracture permeability tensor. The mass conservation equation reads

$$ \nabla_S \cdot j = 0, \qquad (19) $$

where $\nabla_S$ is the 2-D gradient operator in the mean fracture plane.

Because of the classical Poiseuille law, the permeability $\sigma$ of a fracture is expected to be of the order of

$$ \sigma \sim \sigma_0 = \frac{b_0^3}{12}. \qquad (20) $$

The dimensionless fracture permeability $\sigma'$ is defined by

$$ \sigma' = \frac{1}{\sigma_0} \sigma. \qquad (21) $$

Note that if electrical conductivity is addressed, the shape of (18) and (19) remains the same: p should be replaced by the electric potential $\psi$; $\sigma$ would be the fracture conductivity tensor.

It is known that natural fractures are more complex than plane channels. However, the flow behaviour due to the complexity of fracture geometry is not the focus of the present study. This topic was specifically addressed by Mourzenko et al. (1995, 2001). In most of this paper, $\sigma$ is taken to be uniform over each fracture and generally identical for all the fractures. A more complex estimation of the fracture permeability will be given in Section 4.2: in addition some calculations with non-constant fracture permeabilities are described and discussed.

These equations must be supplemented with non-flux conditions at the fracture edges and conservation (for the flux) and continuity (for pressure) equations along the fracture intersections.

Any standard overall boundary condition can be applied to the network. For instance, pressures or fluxes could be applied along inlet and outlet lines drawn on fractures of the network. In the case

where the fracture network can be considered as statistically homogeneous at the field scale, which is assumed to be large with respect to the lateral dimensions R of the fractures, a macroscopic pressure gradient $\overline{\nabla p}$ induces an average flux $\overline{v}_n$, which is related to the pressure gradient by Darcy's law (Adler 1992)

$$ \overline{v}_n = -\frac{1}{\mu} \overline{\nabla \cdot \overline{K}} \cdot \overline{\nabla p}. \qquad (22) $$

$\overline{K}$ is the network permeability tensor [$L^2$]. For isotropic networks, it is a spherical tensor

$$ \overline{K} = \overline{K} \cdot I. \qquad (23) $$

Since the generated numerical samples are spatially periodic in the present case, it is a simple matter to impose a macroscopic pressure gradient $\overline{\nabla p}$ on the infinite periodic medium, and to derive the corresponding components of $\overline{K}$ from the mean flux $\overline{v}_n$.

It is convenient to introduce a dimensionless permeability tensor $\overline{K}'$,

$$ \overline{K}' = \frac{1}{K_0} \overline{K}, \qquad K_0 = \frac{\sigma_0}{R}. \qquad (24) $$

Koudina et al. (1998) developed general numerical tools to solve the local flow equations in fracture networks described by a 3-D triangular mesh, and systematically investigated the flow properties of the same class of random networks of plane fractures as Huseby et al. (1997). They showed that the permeability varies as

$$ \overline{K}' = 0.0455(\rho' - \rho'_c)^{1.57} \quad (3.5 \le \rho' \le 20). \qquad (25) $$

Snow (1969) considered networks where all the fractures are infinite plane channels with an arbitrary orientation distribution. This is equivalent to assuming that the whole surface of all the fractures in the network may contribute to the flow and can be valid only in the limit of very dense networks. For an isotropic network, the permeability tensor is given by

$$ K_{Sn}^{iso} = \frac{2}{3} \sigma S I, \qquad (26) $$

where S is the volumetric surface area of fractures, i.e. the inverse of the length Lf in eq. (6). This result is easily generalized for anisotropic networks by introducing the fracture orientation distribution, which yields a non-spherical tensor $K_{Sn}$

$$ K_{Sn} = \int \int_{\Omega} \sigma(n)S(n)(I - nn) d^2n, \qquad (27) $$

where $\Omega$ is the unit sphere and $S(n)d^2n$ is the volumetric surface area of fractures with normal vector in the solid angle $d^2n$ around n, with permeability $\sigma(n)$.

For finite polygons, S can be expressed in terms of the surface A and perimeter P of the polygons. Hence, the dimensionless permeability in Snow's (1969) model network with the same surface density is

$$ K_{Sn}^{iso'} = \frac{4}{3} \frac{R}{P} \rho', \qquad (28) $$

where R/P is a shape factor, equal, for instance, to $\sqrt{2}/8$, 1/6, and 1/2$\pi$ for square, hexagonal, and circular fractures, respectively. The numerical calculations of Koudina et al. (1998) showed that (28) is indeed a possible asymptote for networks of finite fractures with very large densities.

This analytical approach can be applied to line data. If the $n_{f,p}$ fractures representing each event f have identical size and shapes,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Geological setting of the Baget area. A north–south cross-section of the surveyed area is displayed in the upper left. A geological map is also displayed where the geological formations north of the Alas fault are not detailed.</figcaption></figure>

and if the volumetric area of each event type is simply $1/L_f$, eq. (27) can be rewritten as a sum over the events:

$$K_{Sn} = \sigma \sum_{f} \frac{n_{f,p}}{L_f} (I - n_f n_f). \qquad (29)$$

It is remarkable that except for $\sigma$, all the quantities in (29) are again geometrical parameters which are readily available from the line surveys. In particular, no assumption regarding the size or shape of the fracture is required since they are assumed to be infinite. As a corollary, eq. (29) yields the same prediction of the permeability of monodisperse and polydisperse networks.

This presentation of the flow properties can be summarized as follows. The permeability of a given network can be determined in three different ways: it can be numerically calculated by the tool developed by Koudina et al. (1998); more conveniently, it can be estimated by (25) or (29).

# 3 FIELD DATA

## 3.1 Geological setting

The Baget watershed (see Fig. 2) is located inside the north Pyrenean zone that was highly deformed during late Cenomanian to Tertiary Pyrenean orogeny induced by the transpressional strike-slip motion of the Iberic and European plates along the north Pyrenean fault (see Choukroune 1992, for a review paper on the tectonics of the Pyrenees). Our study concerns the karstified part of the basin consisting of a slice of alternating metamorphic Jurassic to Cretaceous dolomites, limestones and calcareous marls, dipping 70° to 90° southwards under the slaty Albian–Cenomanian Ballongue flysch, remnant of a Cretaceous pull-apart basin opened during strike-slip motion along the Pyrenean margin (Johnson & Clarence 1989).

The Baget drainage basin is limited to the north by the Alas vertical fault, a satellite of the north Pyrenean fault, running mainly west–east. The original stratification is easily seen at the outcrop as ubiquitous open discontinuities running mainly east–west. The same direction is also recognized as the cleavage direction in Ballongue flysch and so produces a major source of anisotropy. A second discontinuity direction is recognized from satellites images, running 170°N to 10°N (Debroas 1987). The present state of stress of the basin is poorly constrained. Goula et al. (1999) proposed that the regional main (compressive) stress is N–S in the eastern Pyrenees from an inversion of focal mechanisms and striation data. However, the world stress map of Reinecker et al. (2003) exhibits only sparse and highly dispersed data around the Baget area. Souriau & Granet (1995) and Souriau et al. (2001) proposed that a rigid block centered around St Gaudens and bordering westward the Baget drainage basin, could explain the lack of seismicity in this area. Hence, the local stress field around the Baget area can only be deduced from local field work. Due to metamorphism, matrix porosity is reduced to less than 1 per cent (Mangin 1974) and voids consist in dissolution caves and in an open fractures and joints. Several caves have been recognized on both sides of the Baget valley. Two of them have been mapped: La Péreyre to the north and St Catherine to the south of the valley. However, the underground drainage system of the Baget valley is far from being completely explored.

Finally, it is worth noting that the river is oriented east–west, a direction which is likely to be followed by most underground waters.

## 3.2 Line surveys

Open cracks orientation and extension have been mapped along three roads: these line surveys correspond to sub-horizontal profiles,
referred to as $P_1$, $P_2$ and $P_3$, oriented 110°, 70° and 120° clockwise from the north, with lengths 700 m, 200 m and 100 m, respectively. In addition, fractured zones consisting in series of subparallel cracks have been characterized by their width and their mean crack spacing.

Some data collected in the La Péreyre and St Catherine cave system have not been included in this survey, but they were used to check the consistency of the data set.

The characteristics of all the recorded events are summarized in Table 1. The x- and y-axes of coordinates are set in the S–N and E–W directions, respectively, and z is vertical and oriented upward. Recall that the vertical component of the unit vectors **n** is taken as positive.

The orientational distribution of the events is shown in Fig. 3. The orientations were measured with a 5° resolution, but they were slightly randomized in Fig. 3 in order to distinguish events with identical orientations. Only very limited information regarding the event extensions can be gained from the present line surveys, since the observed trace lengths are very often truncated by the boundaries of the outcrop along the roads as shown in Table 1.

Fig. 3 clearly shows that the events can be categorized into four main families, which are referred to hereafter as $F_1$ to $F_4$ (see Table 1). Two of them ($F_1$ and $F_4$) are sub-vertical, and roughly in the E–W and N–S directions, whereas $F_2$ and $F_3$ have a slope of about 40°–50°. In addition, a few events do not belong to any of the four families. They are denoted by zeros in the last column of Table 1.

<table><caption>**Table 1.** Fracturation event characteristics from the line surveys: *P* and *x* are the profile number and the location along the profile; *L*<sub>*f*</sub> is the total corrected profile length according to formula (31); *n* is the normal vector to the event plane; *w* is the event width along the profile and *d* is the typical spacing of fractures within the event (*d* = 0 denotes an isolated fracture); *L*<sub>*t*</sub> is the apparent trace length (+ indicates that the traces were truncated by the outcrop boundaries); and *F* is the family that the event belongs to. All distances are in metres.</caption><thead><tr><th></th><th><i>P</i></th><th><i>x</i></th><th><i>L</i><sub>f</sub></th><th><b><i>n</i></b></th><th><i>w</i></th><th><i>d</i></th><th><i>L</i><sub>t</sub></th><th><i>F</i></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>5</td><td>806.31</td><td>-0.1632</td><td>0.3420</td><td>10.00</td><td>0.50</td><td>1.00</td><td>1</td></tr><tr><td>2</td><td>1</td><td>85</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>0.50</td><td>2.00 +</td><td>1</td></tr><tr><td>3</td><td>1</td><td>145</td><td>838.36</td><td>0.0000</td><td>1.0000</td><td>0.0000</td><td>0.20</td><td>2.00 +</td><td>1</td></tr><tr><td>4</td><td>1</td><td>155</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>7.00</td><td>0.20</td><td>2.00 +</td><td>1</td></tr><tr><td>5</td><td>1</td><td>180</td><td>838.36</td><td>0.0000</td><td>1.0000</td><td>0.0000</td><td>20.00</td><td>1.00</td><td>2.00 +</td><td>1</td></tr><tr><td>6</td><td>1</td><td>200</td><td>526.37</td><td>0.0560</td><td>-0.6403</td><td>0.7660</td><td>2.00</td><td>0.50</td><td>2.00 +</td><td>2</td></tr><tr><td>7</td><td>1</td><td>205</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>3.00</td><td>0.50</td><td>2.00 +</td><td>1</td></tr><tr><td>8</td><td>1</td><td>218</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>1.00</td><td>0.50</td><td>2.00 +</td><td>1</td></tr><tr><td>9</td><td>1</td><td>230</td><td>254.40</td><td>0.4924</td><td>-0.4132</td><td>0.7660</td><td>3.00</td><td>1.00</td><td>2.00 +</td><td>3</td></tr><tr><td>10</td><td>1</td><td>245</td><td>838.36</td><td>0.0000</td><td>-1.0000</td><td>0.0000</td><td>50.00</td><td>2.00</td><td>2.00 +</td><td>1</td></tr><tr><td>11</td><td>1</td><td>290</td><td>93.80</td><td>0.6040</td><td>-0.2198</td><td>0.7660</td><td>10.00</td><td>~0.10-1.00</td><td>~0.10-2.0</td><td>3</td></tr><tr><td>12</td><td>1</td><td>295</td><td>838.36</td><td>0.0000</td><td>-1.0000</td><td>0.0000</td><td>50.00</td><td>1.00</td><td>2.00</td><td>1</td></tr><tr><td>13</td><td>1</td><td>300</td><td>489.44</td><td>0.1664</td><td>-0.6209</td><td>0.7660</td><td>5.00</td><td>1.00</td><td>2.00 +</td><td>2</td></tr><tr><td>14</td><td>1</td><td>320</td><td>509.85</td><td>0.1116</td><td>-0.6330</td><td>0.7660</td><td>5.00</td><td>1.00</td><td>2.00</td><td>2</td></tr><tr><td>15</td><td>1</td><td>335</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>30.00</td><td>~1.00-2.00</td><td>1.00</td><td>1</td></tr><tr><td>16</td><td>1</td><td>345</td><td>316.32</td><td>-0.2500</td><td>0.4330</td><td>0.8660</td><td>20.00</td><td>~0.10-2.00</td><td>5.00</td><td>0</td></tr><tr><td>17</td><td>1</td><td>350</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>10.00</td><td>1.00</td><td>2.00 +</td><td>1</td></tr><tr><td>18</td><td>1</td><td>370</td><td>627.30</td><td>0.0668</td><td>-0.7631</td><td>0.6428</td><td>10.00</td><td>1.00</td><td>2.00</td><td>2</td></tr><tr><td>19</td><td>1</td><td>385</td><td>521.57</td><td>0.3237</td><td>-0.6943</td><td>0.6428</td><td>10.00</td><td>2.00</td><td>2.00</td><td>2</td></tr><tr><td>20</td><td>1</td><td>450</td><td>484.63</td><td>0.3830</td><td>-0.6634</td><td>0.6428</td><td>10.00</td><td>0.50</td><td>2.00</td><td>2</td></tr><tr><td>21</td><td>1</td><td>470</td><td>314.10</td><td>-0.8627</td><td>-0.0755</td><td>0.5000</td><td>138.00</td><td>0.50</td><td>2.00 +</td><td>0</td></tr><tr><td>22</td><td>1</td><td>515</td><td>128.70</td><td>0.4330</td><td>-0.2500</td><td>0.8660</td><td>0.00</td><td>0.00</td><td>10.00</td><td>3</td></tr><tr><td>23</td><td>1</td><td>545</td><td>858.14</td><td>0.2588</td><td>0.9659</td><td>0.0000</td><td>2.00</td><td>0.50</td><td>2.00 +</td><td>1</td></tr><tr><td>24</td><td>1</td><td>555</td><td>858.14</td><td>0.2588</td><td>0.9659</td><td>0.0000</td><td>2.00</td><td>0.50</td><td>2.00 +</td><td>1</td></tr><tr><td>25</td><td>1</td><td>570</td><td>128.70</td><td>0.4330</td><td>-0.2500</td><td>0.8660</td><td>0.00</td><td>0.00</td><td>60.00</td><td>3</td></tr><tr><td>26</td><td>1</td><td>590</td><td>639.80</td><td>0.3971</td><td>-0.8517</td><td>0.3420</td><td>10.00</td><td>1.00</td><td>2.00</td><td>0</td></tr><tr><td>27</td><td>1</td><td>600</td><td>838.36</td><td>0.0000</td><td>1.0000</td><td>0.0000</td><td>10.00</td><td>0.10</td><td>1.00</td><td>1</td></tr><tr><td>28</td><td>1</td><td>625</td><td>544.64</td><td>0.5390</td><td>-0.7698</td><td>0.3420</td><td>50.00</td><td>2.00</td><td>4.00 +</td><td>0</td></tr><tr><td>29</td><td>1</td><td>675</td><td>544.64</td><td>0.5390</td><td>-0.7698</td><td>0.3420</td><td>10.00</td><td>1.00</td><td>1.00</td><td>0</td></tr><tr><td>30</td><td>2</td><td>15</td><td>838.36</td><td>0.0000</td><td>1.0000</td><td>0.0000</td><td>5.00</td><td>0.05</td><td>5.00 +</td><td>1</td></tr><tr><td>31</td><td>2</td><td>45</td><td>858.14</td><td>0.2588</td><td>0.9659</td><td>0.0000</td><td>30.00</td><td>0.50</td><td>5.00 +</td><td>1</td></tr><tr><td>32</td><td>2</td><td>70</td><td>409.72</td><td>-0.2500</td><td>-0.4330</td><td>0.8660</td><td>2.00</td><td>0.50</td><td>2.00</td><td>0</td></tr><tr><td>33</td><td>2</td><td>75</td><td>462.28</td><td>0.9397</td><td>0.3420</td><td>0.0000</td><td>0.00</td><td>0.00</td><td>10.00</td><td>4</td></tr><tr><td>34</td><td>2</td><td>80</td><td>323.62</td><td>1.0000</td><td>0.0000</td><td>0.0000</td><td>40.00</td><td>0.50</td><td>5.00 +</td><td>4</td></tr><tr><td>35</td><td>2</td><td>110</td><td>323.62</td><td>1.0000</td><td>0.0000</td><td>0.0000</td><td>10.00</td><td>0.50</td><td>5.00 +</td><td>4</td></tr><tr><td>36</td><td>2</td><td>115</td><td>462.28</td><td>0.9397</td><td>0.3420</td><td>0.0000</td><td>0.00</td><td>0.00</td><td>10.00</td><td>4</td></tr><tr><td>37</td><td>2</td><td>132</td><td>509.85</td><td>0.1116</td><td>-0.6330</td><td>0.7660</td><td>3.00</td><td>1.00</td><td>2.00 +</td><td>2</td></tr><tr><td>38</td><td>2</td><td>155</td><td>462.28</td><td>0.9397</td><td>0.3420</td><td>0.0000</td><td>0.00</td><td>0.00</td><td>5.00</td><td>4</td></tr><tr><td>39</td><td>2</td><td>185</td><td>362.86</td><td>0.2868</td><td>-0.4967</td><td>0.8192</td><td>5.00</td><td>1.00</td><td>10.00 +</td><td>2</td></tr><tr><td>40</td><td>3</td><td>35</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>30.00</td><td>0.10</td><td>2.00 +</td><td>1</td></tr><tr><td>41</td><td>3</td><td>50</td><td>396.59</td><td>-0.0868</td><td>0.4924</td><td>0.8660</td><td>2.00</td><td>1.00</td><td>1.00</td><td>0</td></tr><tr><td>42</td><td>3</td><td>60</td><td>323.62</td><td>1.0000</td><td>0.0000</td><td>0.0000</td><td>20.00</td><td>0.10</td><td>1.00 +</td><td>4</td></tr><tr><td>43</td><td>3</td><td>60</td><td>858.06</td><td>0.1736</td><td>0.9848</td><td>0.0000</td><td>20.00</td><td>0.10</td><td>1.00 +</td><td>1</td></tr><tr><td>44</td><td>3</td><td>60</td><td>419.18</td><td>0.0000</td><td>-0.5000</td><td>0.8660</td><td>20.00</td><td>1.00</td><td>2.00 +</td><td>2</td></tr><tr><td>45</td><td>3</td><td>100</td><td>323.62</td><td>1.0000</td><td>0.0000</td><td>0.0000</td><td>10.00</td><td>0.10</td><td>2.00 +</td><td>4</td></tr><tr><td>46</td><td>3</td><td>100</td><td>396.59</td><td>0.0868</td><td>-0.4924</td><td>0.8660</td><td>0.00</td><td>0.00</td><td>10.00</td><td>2</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Orientation of the fracturation events in Table 1. Symbols correspond to the position of the vector **n** on the unit sphere. Families 1, 2, 3 and 4 are denoted by Δ, ◦, ♦ and ∇, respectively. The line spacing is 10° for the radial angle and 5° for the inclination.</figcaption></figure>

F1 results from the original stratification and from the strike- slip motion that characterizes the Pyrenean orogeny. A direction near 10°N close to F4 has been observed during unpublished microtectonic works of Mangin and Paredes in the Baget basin. In addition, F1, F2 and F4 are observed in the La Péreyre and St Catherine caves together with numerous cracks of different orientations that are not observed at the surface, and which could be related to the mechanical heterogeneity induced by the cave itself.

## 3.3 Data analysis

### 3.3.1 Spatial distribution

For the sake of the subsequent simulations, it is important to determine whether the fracturation events are spatially correlated. Although the data set is too limited for a detailed analysis, some information can be obtained from the 13 events of family $F_1$ on profile $P_1$. The spacings $s_i$ between successive events $i$ and $i + 1$ have an average $\langle s \rangle = 43.3$ m and a standard deviation $\sigma_s = 48.6$ m. Recall that if the events are Poissonian, i.e. without any spatial correlation (see Kingman 1993), the spacing between their intersections with a scan line obeys an exponential probability law proportional to $\exp(-s/\langle s \rangle)$, with $\langle s \rangle = \sigma_s$ (Adler & Thovert 1999). In the present case, $\sigma_s$ is slightly larger than the mean spacing, which indicates that the events are slightly more clustered than in a random distribution.

The variogram of the spacing can be defined as the average (Long & Billaux 1987)

$$ \gamma_s(n) = \langle (s_i - s_{i+n})^2 \rangle. \qquad (30) $$

For a Poisson distribution, the spacings are uncorrelated and $\gamma_s$ is constant and equal to $\langle s \rangle^2$. The variogram $\gamma_s$ is plotted in Fig. 4 for family $F_1$ in profile $P_1$. Deviations from $\langle s \rangle^2$ are observed. However, they are due in most part to the small size of the statistical data set. For instance, if the spacing between the last observed event and the end of the profile is included in the calculation, the deviations are significantly reduced. In addition, this calculation does not take into account the very different widths of the various events.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Variogram (30) of the spacings of the intersections of events from family $F_1$ with profile $P_1$ without (—○—) or with (···○··) inclusion of the part of the profile following the last observed intersection. The thick solid line corresponds to uncorrelated event locations.</figcaption></figure>

No definite intercorrelation was found either between the event width and location, or spacing.

Therefore, in the absence of any conclusive evidence of spatial correlations, it is assumed in the following that the event location and other characteristics are uncorrelated.

### 3.3.2 Volumetric density

The volumetric density of each of the 46 events in Table 1 is given by (6), where $A_f$ ($f = 1$ to 46) is the event area and $L_f$ is given by an obvious extension of (6) in order to take into account the measurements along the three profiles

$$ L_f = \sum_{j=1,3} L_j |\mathbf{p}_j \cdot \mathbf{n}_f|, \qquad (31) $$

where $L_j$ and $\mathbf{p}_j$ ($j = 1, 2, 3$) are the length and orientation of profile $P_j$, respectively. The length $L_f$ is the reduced total length of profile in the direction normal to each event and is given in Table 1.

It is worth noting that although families $F_1$ and $F_2$ yield the most numerous events in the data set, this is mainly due to their favourable orientations relative to the profile axes. In the directions normal to their mean planes, the four families have similar frequencies, with 21, 21, 26 and 18 intersections per kilometre of profile for $F_1$ to $F_4$, respectively.

Some information can be given on the area $A_f$. In addition to the event extensions compiled in Table 1, complementary observations in a quarry in the same area showed that the traces of the fractures generally extend throughout a 25 m high subvertical wall. Moreover, it is believed that joints along the E–W direction, i.e. parallel to the main regional tectonic feature in the Pyrénées, have a larger extension than along the S–N direction. These indications will be used as guidelines in the simulations.

## 3.4 Numerical reconstruction

Numerical samples of fractured rock are stochastically generated according to the two-step procedure described in Section 2.3. The cell size is generally equal to $L = 600$ m.
The events and the fractures they contain are assumed to have identical plane circular or polygonal shapes circumscribed by circles of radius R. In the absence of precise experimental information, two simple cases only are considered here: in the first case, all events have the same size R; in the second case, the events are oriented E-W with the exception of events 34 and 35 that have a size R<sub>M</sub> twice as large as the R<sub>m</sub> of the other events. The larger events are essentially those in family F<sub>4</sub>, and events 11 and 21. Hereafter, these two cases are referred to as monodisperse and bidisperse, respectively.

It is important to note that the assumption about the size in the second case is in agreement with large scale geological observation, since only F1 and F4 are observed on satellite images or on regional geological maps.

Given the size and shape of an event *f*, and thus its area *A*<sub>*f*</sub>, the volumetric density ρ<sub>*f*</sub> results from (31) and the number *n*<sub>*f*,p</sub> to be inserted in the unit cell is ρ<sub>*f*</sub>L<sup>3</sup>. A number of events equal to the integer part of *n*<sub>*f*</sub> is inserted first, and then an additional one with a probability equal to the fractional part of *n*<sub>*f*</sub>. According to the discussion in paragraph 3.3.1, they are placed with their centers at random locations, without any spatial correlation.

Various values of R have been considered, from 25 m to 100 m.

In the second step, each event is replaced by series of *n*<sub>*f*,p</sub> parallel-lel fractures, with the same size and shape, and their centers evenly distributed on the normal axis of the event. *n*<sub>*f*,p</sub> can be derived from the width *w* and fracture spacing *d* of the event, according to (1). Since this can yield up to 300 fractures per event, *f* = 40, *n*<sub>*f*,p</sub>, for instance for, was rescaled by a constant factor so that it does not exceed a maximal value 1 + *n*<sub>max</sub>. This is done by introducing a factor β, equal for all the events, in the definition of *n*<sub>*f*,p</sub>:

$$n_{f,p} = 1 + \text{Int}\left(\frac{w_f}{\beta d_f}\right), \quad \beta = \frac{300}{n_{\text{max}}}. \qquad (32)$$

In order to check that the numerical samples actually comply with the field observations, 500 random realizations were generated, and line surveys were simulated through them. The scan line is 500 m long, oriented 120° clockwise from north, in a horizontal plane at z<sub>0</sub> = L/4 in the samples. The intersections with each type of event were recorded and averaged over the 500 realizations.

Monodisperse networks of circular events were considered first. Examples of trace maps in a horizontal section are shown in Fig. 5, with the events represented either by single fractures or by series of parallel fractures. The number *n*<sub>*f,i*</sub> of event intersections with the scan line was checked in the first type of network. It is compared to the value expected from the scan line length, orientation and from the event parameters *L*<sub>*f*</sub> and *n*<sub>*f*</sub> in Fig. 6. The agreement is quite satisfactory. Note that the statistical error bar is of the order of $\sqrt{n_{f,i}(1 - n_{f,i})/500} \sim 0.02$.

Fig. 7 shows examples of bidisperse networks. In this case, the event radius is R<sub>M</sub> = 100 m for those oriented in the W-E direction, and R<sub>m</sub> = 50 m for the others. All other parameters are identical to the previous case, except for the number of events, since it depends on their size (see eq. 6). The comparison of the mean numbers of event intersections with the scan line in the simulated networks with the expected values was as good as in the monodisperse case.

For illustration, Fig. 8 shows a 3-D view of a bidispersed reconstructed sample, with R<sub>M</sub> = 50 m for the W-E events, R<sub>m</sub> = 25 m for the others, and L = 200 m. The events are represented by single hexagonal fractures.

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

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Trace maps obtained in an horizontal section through monodisperse reconstructed samples. The cell size is L = 600 m and the event radius R = 100 m. The events are represented by single fractures in (a) and as series of fractures in (b). The solid lines correspond to the traces of the events which intersect the 120°N, 500 m long scan line (thick solid line). Traces in dotted lines do not intersect the scan line. Distances are in *R* units, i.e. 1:100 m.</figcaption></figure>

## 3.5 Connectivity

The connectivity can now be easily quantified by the formulae derived in Section 2.4. The sums in (11) and (12) can be directly calculated from the line survey data, without any hypothesis regarding
<figure><img src="image_9.png" /><figcaption>Figure 6. Mean number of event intersections with a 500 m long, 120°N scan line (see Fig. 5a). The abscissae correspond to the event number given in the first row of Table 1. Data in (a) correspond to the average over 500 monodisperse reconstructed samples, with L = 600 m, R = 50 m. Data in (b) are the expected values resulting from the field data.</figcaption></figure>

the event size and shape. They are equal to

$$ \sum_{i} \sum_{j} \frac{\sin y_{i,j}}{L_{i} L_{j}} = 8.081510^{-3} \, \text{m}^{-2} \quad \sum_{i} \frac{1}{L_{i}} = 0.1078 \, \text{m}^{-1}. \quad (33) $$

The mean number $\rho'$ of intersections per event is thus expressed as

$$ \rho' = \frac{2\rho_I}{\rho} = 0.0477P, \quad (34) $$

where P in the last equality is expressed in metres. This result shows that the assumption made for the event shape in the modelization has little influence on the connectivity, since, with R in metres, (13) yields

$$ \begin{aligned} \rho' &= 0.300R(\text{discs}), \quad \rho' = 0.286R(\text{hexagons}), \\ \rho' &= 0.270R(\text{squares}). \end{aligned} \quad (35) $$

The global effective exclusion volume $V_{ex}$ can be derived from (14)

$$ V_{ex} = \frac{\rho'}{\rho} = 0.4424AP, \quad (36) $$

a value which is slightly smaller than AP/2 obtained for isotropic networks.

Hence, the present anisotropic distribution of the event orientations tends to decrease the connectivity, in terms of mean number of connections per event given the global volumetric density, with respect to the isotropic case. However, this is only a global statement, which does not take into account the fact that events belonging to

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

<figure><img src="image_11.png" /><figcaption>Figure 7. Trace maps obtained in a horizontal section through bidisperse reconstructed samples. The cell size is L = 600 m and the event radius $R_M = 100$ m for those oriented W-E, and $R_m = 50$ m for the others. The events are represented by a single fractures in (a) and as series of fractures in (b). The solid lines correspond to the traces of the events which intersect the 120°N, 500 m long scan line (thick solid line). Traces in dotted lines do not intersect the scan line. Distances are in $R_M$ units, i.e. 1:100 m.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Example of bidisperse reconstructed sample, with $R_M = 50$ m for the W–E events, $R_m = 25$ m for the others, and $L = 200$ m. The network contains 434 events represented by single fractures. Coordinates are in $R_M$ units, i.e. 1:50 m. The density $\rho'$ of the network is equal to 7.60.</figcaption></figure>

the same family are unlikely to intersect, whereas events from fam-
ilies $F_1$ and $F_2$, which are nearly orthogonal, have a much greater
probability of intersection.

**3.6 Probability of percolation**

The numerical tools applied in this Section were originally devised
by Huseby *et al.* (1997), who gave a full description of the theoretical
framework and methods in use. Since, they have been significantly
improved in terms of computational efficiency, but without any ma-
jor change in the main features.

Let us insist again on the fact that the main focus of this paper is the
analysis of the connectivity as it is controlled by the dimensionless
density ρ'; the role of the individual fracture permeability will be
briefly addressed in Section 4.2.

The percolation of the reconstructed samples was investigated
for hexagonal fractures, with circumscribed disc radius R. As pre-
viously mentioned (see eq. 35), the particular choice of the fracture
shape has a limited influence on the network connectivity: more
important is the size of the events, since it directly conditions the
volumetric density ρ<sub>f</sub>, which corresponds to the intersection spacing
(or volumetric area) observed on the field.

The probability of percolation Π<sub>p</sub> of these networks was studied
as a function of the density ρ' and of the relative sample size L/R,
since it is well-known that close to the percolation threshold, many
properties of the medium, including Π<sub>p</sub>, are size-dependent (Fisher
1971); instead of the sudden transition to percolation at the threshold
of infinite systems, the probability of percolation of finite systems
gradually increases from zero to one, as ρ' increases. The transition
becomes sharper as the system size increases.

Percolation was studied along the x (S–N) and y (E–W) directions
assuming periodic boundary conditions in all directions. Recall that
because of the orientation of the Lachein river, flow is likely to be
east–west in the karstified limestone. The relative cell size L/R was
varied from 4 to 16, and the density ρ' was gradually increased. The
variations of ρ' can be interpreted either as variations of the event
size R with the field data kept unchanged (see eq. 13, with P = 6R
for hexagons), or as variations of the density for a constant size, with
all the lengths Lᵢ multiplied by a common coefficient in eq. (13). A
number N = 500 of random realizations were generated for each
value of the parameters, and the probability of percolation Πₚ was
defined as the fraction of percolating configurations.

The critical density ρ'_{LC} is defined as the value for which Π_p(ρ')
is equal to 1/2. It is determined by fitting the data for Π_p(ρ') with
<table><caption>**Table 2.** Parameters of the fitted scaling laws 38.</caption><thead><tr><th>Configuration</th><th>n<sub>max</sub></th><th>Direction</th><th>ρ'<sub>c</sub></th><th>b</th><th>ν</th></tr></thead><tbody><tr><td>Mono</td><td>0</td><td>X, Y, Z</td><td>2.254</td><td>3.8165</td><td>0.854</td></tr><tr><td>Mono</td><td>10</td><td>X, Y, Z</td><td>3.004</td><td>5.5085</td><td>0.940</td></tr><tr><td>Bi</td><td>0</td><td>X</td><td>1.951</td><td>2.7567</td><td>1.068</td></tr><tr><td>Bi</td><td>0</td><td>Y</td><td>1.984</td><td>4.5424</td><td>0.888</td></tr><tr><td>Bi</td><td>0</td><td>Z</td><td>2.026</td><td>5.2961</td><td>0.839</td></tr><tr><td>Bi</td><td>10</td><td>X</td><td>3.014</td><td>9.1854</td><td>0.883</td></tr><tr><td>Bi</td><td>10</td><td>Y</td><td>3.117</td><td>8.5922</td><td>0.932</td></tr><tr><td>Bi</td><td>10</td><td>Z</td><td>3.031</td><td>8.7533</td><td>0.913</td></tr></tbody></table>

an error function of the form (Stauffer & Aharony 1994; Reynolds et al. 1980):

$$ \Pi_{\rho}(\rho') = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\rho'} \frac{1}{\Delta_L} \exp \left\{ -\left[ \frac{\xi - \rho'_{Lc}}{2\Delta_L} \right]^2 \right\} d\xi, \quad (37) $$

where $\Delta_L$ and $\rho'_{Lc}$ are fit parameters. Such a function provides an excellent fit to the data obtained for monodisperse and bidisperse networks with the presence or absence of parallel fractures.

Once $\Delta_L$ and $\rho'_{Lc}$ have been evaluated for several values of L, the asymptotic value of the density $\rho'_c$ for infinite systems can be derived from the two scaling relations (Fischer 1971; Charlaix 1986):

$$ \rho'_{Lc} - \rho'_c \propto L^{-1/\nu} \quad \Delta_L \sim bL^{-1/\nu}. \quad (38) $$

Plots of $\rho'_{Lc}$ as a function of $\Delta_L$ were extrapolated to zero in order to find $\rho'_{c}$. Similarly, the critical exponent $\nu$ was obtained by plotting $\Delta_L$ as a function of $L/R$. The data obtained are summarized in Table 2.

The results for monodisperse samples, with fracturation events represented by a single fracture, are shown in Fig. 9(a), where the sharpening of the transition to percolation as $L/R$ increases is clearly visible. The critical density $\rho'_c$ as $\Delta_L \to 0$ was found to be $\rho'_c = 2.254$, which is close to the result $\rho'_c \approx 2.26$ of Huseby et al. (1997) for randomly-oriented fractures. In view of (13), the fracturation in the Baget basin has a density $\rho' \sim 2.4$ if the event size is $R \sim 8.3$ m. Since trace lengths much longer than this value have been observed many times when the outcrop geometry allowed for it, i.e. in the quarry or when the traces were roughly parallel to the roads, the connectivity of the fracture network is probably very large and it percolates.

It is nevertheless interesting to investigate the percolation in this low-density range, in order to assess the influence of size polydispersity and of the width of the fracturation events.

Fig. 9(b) shows results obtained with bidisperse networks. The percolation was checked in the x (S–N) and y (E–W) directions. The latter is parallel to the largest events and to the river, with a radius $R_M$, whereas the former is in the orthogonal direction. The cell size $L/R_M$ was varied from 4 to 8, causing the same sharpening of the transition to unit percolation of probability as in the monodisperse case.

Recall that the fracturation is not only bidisperse, but also anisotropic, since the largest events are all aligned in the E–W direction. The critical densities $\rho'_{Lc}$ measured along the two axes are indeed slightly different; they are equal to 2.17 (S–N) and 2.10 (E–W), for $L/R_M = 8$. However, this might result, at least partly, from different size effects in the x and y directions. The asymptotic values found as $\Delta_L \to 0$ (Table 2) are indeed very close in all the directions. They were respectively $\rho'_c = 1.951$ and 1.984 in the x and y directions. In any case, both values are smaller than for the monodisperse case. This means that this anisotropic polydispersity improves the efficiency of event connections to ensure global percolation.

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Probability of percolation $P_p$ in the x direction as a function of the density $\rho'$, for monodisperse (a) and bidisperse (b) networks of fractures. The data are averages over 500 random realizations. In the bidisperse networks (b) the size $R_M$ of the east–west fractures is twice the size $R_m$ of the others. The symbols $\diamond$, $\circ$, $\nabla$ and $\square$ correspond to cell sizes $L/R = 16, 12, 10$ and 8, respectively in (a), and to $L/R_M = 8, 6, 5$ and 4 in (b).</figcaption></figure>

The previous calculations were repeated for fracturation events with non-zero width, represented by series of up to $n_{max}$ parallel fractures, according to the observed event width $w_f$.

In this situation, the network connectivity can be quantified by two dimensionless densities, namely the mean number $\rho'_e$ of intersections per event, and the mean number $\rho'_p$ of fracture intersections per fracture. If all events are single fractures, $\rho'_e = \rho'_f$. If $w_f > 0$, $\rho'_e$ can be easily estimated from the field data by modifying (8), as noted in Section 2.4. On the other hand, it is difficult to generalize the argument in Section 2.4 to predict $\rho'_p$, since it only applies for objects without spatial correlations. In the reconstructed samples, $\rho'_p$ is determined according to (13), and $\rho'_c$ is determined likewise, by considering that two events are intersecting if they contain a pair of intersecting fractures.

In addition, two different percolation criteria can be introduced. One may assume that all the fractures within an event are mutually
connected: this may represent a situation where the whole event is a high permeability region. On the other hand, one may consider that the series of parallel fractures are unconnected, because they are separated by layers of intact rock. Then, a cluster of overlapping events may actually not contain a spanning cluster of connected fractures. The latter point of view was adopted here, and the percolation probability is analyzed in terms of the density $\rho'_{p}$. The subscript $p$ is omitted in order to simplify the notations, but it should be remembered that $\rho'_{p}$ stands for $\rho'_{p}$.

Size effects are much stronger than for networks of single fractures: substantial variations of $\rho'_{c}$ still take place between $L/R = 12$ and $L/R = 16$, for the monodisperse case, and between $L/R_{M} = 6$ and $L/R_{M} = 8$ for the bidisperse case. Therefore, the asymptotic value $\rho'_{c}(\Delta_{L} \to 0)$ were estimated for networks containing series of parallel fractures. For a monodisperse configuration, we found $\rho'_{c} = 3.042$, which is a much larger value than the results found when we did not consider parallel fractures. For the bidisperse configuration with parallel fractures, we found $\rho'_{c} = 3.014$ and 3.117 in the x and y directions. These densities are clearly larger than the result $\rho'_{c} \sim 2$ obtained when the events are represented by a single fracture. This is an expected result, since $\rho'$ is evaluated from the total number of fracture intersections. When two series of parallel fractures intersect, many fracture intersections, which are largely redundant for the network percolation, are created. The number of such intersections is roughly a quadratic function of the number of fractures in the series, whereas the total number of fractures is a linear function. Thus, $\rho'$ should increase roughly linearly with the size of the series.

This is confirmed by the data in Fig. 10, which summarizes the results for large cell size $L/R = 16$ for monodisperse networks and $L/R_{M} = 8$ for bidisperse networks. The fracturation events are represented by a single fracture, or by series of up to $n_{max} = 10$ or 20 parallel fractures. The corresponding curves for $\Pi_{p}$ are shifted towards larger densities $\rho'$, i.e. larger numbers of intersections per fracture, when the size of the series increases; furthermore, the shift for $n_{max} = 20$ is about twice the shift for $n_{max} = 10$.

## 3.7 Flow properties

The flow equations described in Section 2.6 were solved in the networks listed in Table 3. The unit cell size $L$ was varied from 200 to 600 m and the event size $R$ from 25 to 100 m. In addition, monodisperse and bidisperse size distributions and various values of $n_{max}$ were considered. As usual, in the bidisperse cases, the size $R_{M}$ of the events in the E–W direction is twice as large as the size $R_{m}$ of the others. Finally, the dimensionless fracture permeability $\sigma'$ was generally taken constant, equal to unity, but in some cases, it was set equal to unity for the events in the E–W direction and to 1/2 for the others. This last case is discussed in Section 4.

As in the previous Section, it was considered that parallel fractures belonging to the same event are mutually unconnected. Thus, when $n_{max} > 1$, the network permeability is analysed in terms of the density of fractures, and $\rho'$ stands for $\rho'_{p}$.

As shown in Table 3, these samples generally contain $ca N_{f} = 500–1000$ of fractures. A 3-D view of one of these samples is displayed in Fig. 8.

The permeability can be evaluated first according to the model elaborated by Snow (1969) (see eq. 27). It is known that this model is valid only in the limit of very dense networks (Koudina et al. 1998), but in view of the dimensionless densities $\rho'$ in Table 3, it can be expected to apply reasonably in the present case. This property is due to the fact that the number of intersections per fracture $\rho'$ increases linearly with the fracture density; when $\rho'$ is large, each fracture

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Probability of percolation $P_{p}$ of monodisperse (______) and bidisperse (---) networks, in the x (S–N, a) and y (E–W, b) directions, as a function of the density $\rho'$. The data are averages over 500 random realizations. In bidisperse networks, the size $R_{M}$ of the east–west fractures is twice the size $R_{m}$ of the others. The cell size is $L/R = 16$ (monodisperse) or $L/R_{M} = 8$ (bidisperse). The events are represented by series of up to $n_{max} = 10$ (○) or 20 (□) parallel fractures, or by single fractures (▽).</figcaption></figure>

behaves as if it were of infinite extension which is precisely the conditions for the application of the Snow formula.

It should be recalled that (29) yields the same result for the permeability of monodisperse and bidisperse networks. For instance, the events that are twice as small in the bidisperse model are four times more numerous (see eq. 6), in order to comply with the observed value of $L_{f}$, and their volumetric area is unchanged.

In view of eq. (32), $K_{Sh}$ can be written as

$$K_{Sh} = \sigma \sum_{f} \frac{1}{L_{f}} \left[ 1 + \text{Int} \left( \frac{n_{\text{max}} w_{f}}{300 d_{f}} \right) \right] (I - n_{f} n_{f}), \quad (39a)$$

$$\approx \sigma \left[ \sum_{f} \frac{1}{L_f} (I - n_f n_f) + n_{\text{max}} \sum_{f} \frac{1}{L_f} \text{Int} \left( \frac{w_f}{300 d_f} \right) (I - n_f n_f) \right]. \quad (39b)$$

The two sums in (39b) are equal to

$$\sum_{f} \frac{1}{L_f} (I - n_f n_f) = 0.0719 \begin{pmatrix} 1.064 & 0.046 & -0.178 \\ 0.046 & 0.996 & 0.230 \\ -0.178 & 0.230 & 0.939 \end{pmatrix} \text{m}^{-1}, \quad (40a)$$
<table><caption>**Table 3.** Permeability calculations: (a) corresponds to monodisperse samples; the effects of the bidispersivity are considered in (b); and the effect of multiple parallel fractures in the events are studied in (c). The cell and event sizes, *L* and *R* or *R*<sub>m</sub>, are given first, in metres. Column *Cfg* describes the configurations, with monodisperse (M) or bidisperse (B) size distributions, and fracture permeabilities either constant (U) or set according to the fracture orientation (V). *n*<sub>max</sub> is the maximum number of parallel fractures in the events. *N*<sub>f</sub>, *ρ* and *ρ'* are the number of fractures, the density in fractures per *R*<sup>3</sup> volume, and the dimensionless density, respectively. *K* is the dimensional mean permeability, with *K*/*σ* in m<sup>-1</sup>, and *K*′ is the dimensionless mean permeability. *k*<sub>i</sub> is the normalized diagonal component of the permeability tensor, *k*<sub>i</sub> = *K*<sub>ii</sub>/K. Subscript *Sn* refers to Snow's model (27).</caption><thead><tr><th></th><th><i>L</i></th><th><i>R</i><sub>M</sub></th><th>Cfg</th><th><i>n</i><sub>max</sub></th><th><i>N</i><sub>f</sub></th><th><i>ρ</i></th><th><i>ρ</i>'</th><th><i>k</i><sub>x</sub></th><th><i>k</i><sub>y</sub></th><th><i>k</i><sub>z</sub></th><th><i>10<sup>2</sup>K</i>/<i>σ</i></th><th><i>K</i>'</th><th><i>k</i><sub>Sn,x</sub></th><th><i>k</i><sub>Sn,y</sub></th><th><i>k</i><sub>Sn,z</sub></th><th><i>10<sup>2</sup>K</i><sub>Sn</sub>/<i>σ</i></th><th><i>K</i>'<sub>Sn</sub></th></tr></thead><tbody><tr><td colspan="18" style="text-align: left; font-weight: bold;">(a)</td></tr><tr><td>1</td><td>200</td><td>25</td><td>M,U</td><td>0</td><td>532</td><td>1.04</td><td>7.15</td><td>1.09</td><td>1.03</td><td>0.93</td><td>2.68</td><td>0.67</td><td>1.07</td><td>0.99</td><td>0.94</td><td>7.20</td><td>1.80</td></tr><tr><td>2</td><td>300</td><td>37.5</td><td>M,U</td><td>0</td><td>802</td><td>1.57</td><td>10.9</td><td>1.04</td><td>1.00</td><td>0.95</td><td>3.89</td><td>1.46</td><td>1.07</td><td>1.00</td><td>0.94</td><td>7.23</td><td>2.71</td></tr><tr><td>3</td><td>200</td><td>50</td><td>M,U</td><td>0</td><td>136</td><td>2.13</td><td>14.7</td><td>1.10</td><td>0.94</td><td>0.96</td><td>4.68</td><td>2.34</td><td>1.07</td><td>0.98</td><td>0.94</td><td>7.36</td><td>3.68</td></tr><tr><td>4</td><td>400</td><td>50</td><td>M,U</td><td>0</td><td>1066</td><td>2.08</td><td>14.0</td><td>1.06</td><td>1.01</td><td>0.94</td><td>4.62</td><td>2.31</td><td>1.06</td><td>1.00</td><td>0.94</td><td>7.22</td><td>3.61</td></tr><tr><td>5</td><td>300</td><td>75</td><td>M,U</td><td>0</td><td>201</td><td>3.14</td><td>21.6</td><td>1.09</td><td>0.97</td><td>0.95</td><td>5.59</td><td>4.19</td><td>1.06</td><td>1.00</td><td>0.94</td><td>7.25</td><td>5.44</td></tr><tr><td>6</td><td>600</td><td>100</td><td>M,U</td><td>0</td><td>899</td><td>4.16</td><td>28.2</td><td>1.05</td><td>1.01</td><td>0.94</td><td>5.90</td><td>5.90</td><td>1.06</td><td>1.00</td><td>0.94</td><td>7.21</td><td>7.21</td></tr><tr><td colspan="18" style="text-align: left; font-weight: bold;">(b)</td></tr><tr><td>1</td><td>200</td><td>25</td><td>M,U</td><td>0</td><td>532</td><td>1.04</td><td>7.15</td><td>1.09</td><td>1.03</td><td>0.93</td><td>2.68</td><td>0.67</td><td>1.07</td><td>0.99</td><td>0.94</td><td>7.20</td><td>1.80</td></tr><tr><td>2</td><td>200</td><td>50</td><td>M,U</td><td>0</td><td>136</td><td>2.13</td><td>14.7</td><td>1.10</td><td>0.94</td><td>0.96</td><td>4.68</td><td>2.34</td><td>1.07</td><td>0.98</td><td>0.94</td><td>7.36</td><td>3.68</td></tr><tr><td>3</td><td>200</td><td>50</td><td>B,U</td><td>0</td><td>434</td><td>6.79</td><td>7.60</td><td>0.99</td><td>1.07</td><td>0.95</td><td>3.28</td><td>1.64</td><td>1.07</td><td>0.99</td><td>0.94</td><td>7.18</td><td>3.59</td></tr><tr><td>4</td><td>200</td><td>50</td><td>B,V</td><td>0</td><td>434</td><td>6.78</td><td>7.60</td><td>0.91</td><td>1.15</td><td>0.96</td><td>2.12</td><td>1.06</td><td>0.96</td><td>1.08</td><td>0.96</td><td>4.46</td><td>2.23</td></tr><tr><td>5</td><td>300</td><td>37.5</td><td>M,U</td><td>0</td><td>802</td><td>1.57</td><td>10.9</td><td>1.04</td><td>1.00</td><td>0.95</td><td>3.89</td><td>1.46</td><td>1.07</td><td>1.00</td><td>0.94</td><td>7.23</td><td>2.71</td></tr><tr><td>6</td><td>300</td><td>75</td><td>M,U</td><td>0</td><td>201</td><td>3.14</td><td>21.6</td><td>1.09</td><td>0.97</td><td>0.95</td><td>5.59</td><td>4.19</td><td>1.07</td><td>1.00</td><td>0.94</td><td>7.25</td><td>5.44</td></tr><tr><td>7</td><td>300</td><td>75</td><td>B,U</td><td>0</td><td>655</td><td>10.2</td><td>11.5</td><td>1.01</td><td>1.04</td><td>0.94</td><td>4.48</td><td>3.36</td><td>1.07</td><td>1.00</td><td>0.94</td><td>7.24</td><td>5.43</td></tr><tr><td>8</td><td>300</td><td>75</td><td>B,V</td><td>0</td><td>655</td><td>10.2</td><td>11.5</td><td>0.92</td><td>1.11</td><td>0.97</td><td>2.84</td><td>2.13</td><td>0.95</td><td>1.08</td><td>0.97</td><td>4.51</td><td>3.38</td></tr><tr><td colspan="18" style="text-align: left; font-weight: bold;">(c)</td></tr><tr><td>1</td><td>200</td><td>25</td><td>M,U</td><td>0</td><td>532</td><td>1.04</td><td>7.15</td><td>1.09</td><td>1.03</td><td>0.93</td><td>2.68</td><td>0.67</td><td>1.07</td><td>0.99</td><td>0.94</td><td>7.20</td><td>1.80</td></tr><tr><td>2</td><td>200</td><td>25</td><td>M,U</td><td>5</td><td>768</td><td>1.50</td><td>10.3</td><td>0.94</td><td>1.08</td><td>0.98</td><td>4.80</td><td>1.20</td><td>0.96</td><td>1.02</td><td>1.03</td><td>10.4</td><td>2.60</td></tr><tr><td>3</td><td>200</td><td>25</td><td>M,U</td><td>10</td><td>1065</td><td>2.08</td><td>14.0</td><td>0.91</td><td>1.08</td><td>1.01</td><td>6.96</td><td>1.74</td><td>0.89</td><td>1.01</td><td>1.09</td><td>14.4</td><td>3.60</td></tr><tr><td>4</td><td>200</td><td>50</td><td>M,U</td><td>0</td><td>136</td><td>2.13</td><td>14.7</td><td>1.10</td><td>0.94</td><td>0.96</td><td>4.68</td><td>2.34</td><td>1.07</td><td>0.98</td><td>0.94</td><td>7.36</td><td>3.68</td></tr><tr><td>5</td><td>200</td><td>50</td><td>M,U</td><td>5</td><td>192</td><td>3.00</td><td>20.6</td><td>0.97</td><td>1.04</td><td>0.99</td><td>7.22</td><td>3.61</td><td>0.95</td><td>1.03</td><td>1.02</td><td>10.4</td><td>5.20</td></tr><tr><td>6</td><td>200</td><td>50</td><td>M,U</td><td>10</td><td>264</td><td>4.13</td><td>27.6</td><td>0.92</td><td>1.09</td><td>0.99</td><td>9.66</td><td>4.83</td><td>0.89</td><td>1.03</td><td>1.09</td><td>14.3</td><td>7.14</td></tr><tr><td>7</td><td>200</td><td>50</td><td>B,U</td><td>0</td><td>434</td><td>6.79</td><td>7.60</td><td>0.99</td><td>1.07</td><td>0.95</td><td>3.28</td><td>1.64</td><td>1.07</td><td>0.99</td><td>0.94</td><td>7.18</td><td>3.59</td></tr><tr><td>8</td><td>200</td><td>50</td><td>B,U</td><td>5</td><td>577</td><td>9.02</td><td>11.5</td><td>0.90</td><td>1.08</td><td>1.02</td><td>5.84</td><td>2.92</td><td>0.96</td><td>1.02</td><td>1.03</td><td>10.5</td><td>5.18</td></tr><tr><td>9</td><td>200</td><td>50</td><td>B,V</td><td>0</td><td>434</td><td>6.78</td><td>7.60</td><td>0.91</td><td>1.15</td><td>0.96</td><td>2.12</td><td>1.06</td><td>0.96</td><td>1.08</td><td>0.96</td><td>4.46</td><td>2.23</td></tr><tr><td>10</td><td>200</td><td>50</td><td>B,V</td><td>5</td><td>577</td><td>9.02</td><td>11.5</td><td>0.80</td><td>1.15</td><td>1.05</td><td>3.90</td><td>1.95</td><td>0.81</td><td>1.13</td><td>1.08</td><td>6.88</td><td>3.44</td></tr></tbody></table>

$$ \sum_f \frac{1}{L_f} \text{Int} \left( \frac{w_f}{300 d_f} \right) (I - n_f n_f) = 0.007985 \begin{pmatrix} 0.682 & -0.054 & 0.132 \\ -0.054 & 0.980 & 0.054 \\ 0.132 & 0.054 & 1.338 \end{pmatrix} \text{m}^{-1}. \quad (40b) $$

The mean permeability $\overline{K}_{Sn}$ over the three axes (one third of the trace) is given by

$$ \overline{K}_{Sn} \approx \sigma(0.0719 + 0.007985n_{\max}) \text{m}^2. \quad (41a) $$

The corresponding dimensionless value is

$$ \overline{K}'_{Sn} \approx \frac{R \overline{K}_{Sn}}{\sigma} = R(0.0719 + 0.007985n_{\max}). \quad (41b) $$

However, it should be noted that the insertion of additional fractures in the events, i.e. the increase of $n_{max}$, also modifies the anisotropy of the network permeability tensor, since the two contributions in eq. (40) have different anisotropies.

Again these different formulae make a clear distinction between the geometric properties that are related to measured quantities, and the individual fracture permeabilities that could be obtained by a different set of measurements. Note that variable fracture permeabilities $\sigma$ can easily be accounted for by keeping $\sigma_f$ under the summation in (29).

The permeabilities obtained from (29) are given in Table 3, for all the investigated cases. The data may slightly differ from (39–41), because $K_{Sn}$ was evaluated on the actual random realizations, instead

of from the data in Table 1, but the differences are always small, both in magnitude and anisotropy: slight differences only occur for the sample containing the smallest number of fractures (136), and they are due to statistical fluctuations. Table 3(b) shows that $K_{Sn}$ is identical for monodisperse and bidisperse networks, and the data in Table 3(c) for $n_{max} > 0$ illustrate the increase in mean permeability and the change in the anisotropy induced by the fracture series.

The permeabilities $\overline{K}$ calculated by solving the 3-D flow equations for monodisperse networks are given in Table 3(a). The two calculations performed for $R = 50$ m ($\rho' \ge 14$) in cells with sizes $L = 200$ and 400 m yield similar results. This agrees with the observation of Koudina et al. (1998) that effects of cell size on the permeability of periodic fracture networks vanish when $\rho' \ge 4$ and $L/R \ge 4$. Hence, $(4R)^3$ samples are large enough to obtain a representative value of the network permeability. As expected, $\overline{K}$ is an increasing function of the network density, quantified by $\rho$ or $\rho'$. It appears that the model (29) overestimates the permeability, which is natural since it assumes that the whole fracture area fully contributes to the flow. However, the error is only of the order of 20–25 per cent for very connected networks ($\rho' \ge 20$), which is a fair result, considering that the model (29) does not require the solution of any flow equation. Moreover, the anisotropy of $\overline{K}$ is very close to that of $K_{Sn}$, even for the smallest density.

The various determinations of the mean permeability are plotted in Fig. 11 as functions of the dimensionless network density $\rho'$. For $\rho' > 10$, $\overline{K}$ increases quasi-linearly with the density. The prediction of (29) yields another straight line, parallel to the former, but shifted vertically. Finally, the result (25) of Koudina et al. (1998)
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Dimensionless mean permeability $\bar{K}'$ of monodisperse networks as a function of the density $\rho'$. Data correspond to the numerical solution of the 3-D flow equations (— ○ —), to the prediction (29) (—), and to the result (25) of Koudina *et al.* (1998) for isotropic networks ( · · ◊ · · ). Data are for: monodisperse networks without parallel fractures (— ○ —), monodisperse networks with parallel fractures ($n_{max}$ = 5, +; and $n_{max}$ = 10, □), bidisperse networks without parallel fractures (×).</figcaption></figure>

for isotropic networks is also presented for comparison. It is in good agreement with the present calculations for small and intermediate densities, but it seems to increase at a faster rate for large densities. The numerical results for bidisperse networks are given in Table 3(b). They give rise to several observations.

Consider first the pairs of cases in lines (2, 3) and (6, 7) in the Table 3(b), which have identical size L, radius $R_M$, constant fracture permeability $\sigma$, but are either monodisperse or bidisperse. Since the volumetric areas of each type of event are identical, (29) results in the same permeability tensor $K_{Sh}$. However, the actual calculations yield different permeabilities for the monodisperse and bidisperse cases. $\bar{K}$ is smaller in the bidisperse networks. This is mostly due to a lesser connectivity, with $\rho'$ twice as small as in the monodisperse case. The permeability anisotropy is also different in the monodisperse and bidisperse cases, with an increase of the permeability in the y direction, parallel to the longer fractures, which is not predicted by the model (29).

Consider next the pairs of cases in lines (1, 3) and (5, 7) in Table 3(b), which have identical size L, constant fracture permeability $\sigma$, but different radii $R_M$, so that their dimensionless densities $\rho'$ are roughly equal. Again, the volumetric areas of each type of event are identical and (29) results in identical permeability tensor $K_{Sh}$. However, the mean permeabilities $\bar{K}$ from the solution of the

flow equations are again different in the monodisperse and bidisperse cases. $\bar{K}$ is now larger in the bidisperse networks. The tensor anisotropies are also different, though in a lesser respect than in the previous comparison.

Thus, as expected, Snow's model (29) poorly succeeds in accounting for the effect of size polydispersity.

The numerical results for monodisperse and bidisperse networks where the fracturation events are represented by series of parallel fractures are given in Table 3(c). In all cases, the presence of parallel fractures increases the permeability of the networks.

The two sets of lines (1–3) and (4–6) in Table 3(c) are series of identical cases, except for increasing values of $n_{max}$, equal to 0, 5 and 10. Accordingly, the number $N_f$ increases, by factors of about 3/2 and 2. Since Snow's model (29) supposes perfect connectivity, $\bar{K}_{Sn}$ increases roughly in the same proportions, although the anisotropy of $K_{Sn}$ changes according to (39). However, the fracture connectivity, quantified by the number of intersections per fracture $\rho'$ also increases with $n_{max}$. For the first set (lines 1–3), $\rho'$ increase from 7.15 to 10.3 and 14.0; therefore, the calculated permeabilities $\bar{K}$ increase in a larger proportion than $N_f$, by factors about 1.8 and 2.6 for $n_{max}$ = 5 and 10, respectively. For the second set (lines 4–6), $\rho'$ is always large (14.7–27.6), and its increase has very little effect. Accordingly, $\bar{K}$ increases by the same factor as $N_f$.

Although the addition of parallel fractures in some of the events increases the permeability, the increment is smaller than when the same number of fractures is added by changing the density of events. Compare for instance line 6 in Tables 3(a) and (c). The two samples have nearly equal densities $\rho$ and $\rho'$; however, $\bar{K}$ is 22 per cent larger when the fractures are evenly distributed among all the event types ($n_{max}$ = 0), than when many of them are concentrated in a few types of events ($n_{max}$ = 10). This is illustrated by the plot of $\bar{K}'$ as a function of $\rho'$ in Fig. 11, for $n_{max}$ = 5 and 10, in comparison with the data for $n_{max}$ = 0 in Table 3(a).

The anisotropy of the permeability tensors also varies when $n_{max}$ increases. The x direction (N–S) is the most favourable to the flow for $n_{max}$ = 0, and becomes the least favourable for $n_{max}$ = 10. Conversely, the y direction becomes preferential. This last feature is unexpected, in view of (40b), which shows that the additional fractures are oriented mostly parallel to the z-axis. Accordingly, the anisotropy of $K_{Sn}$ moves from the x to the z direction as $n_{max}$ increases. The different behaviour of $\bar{K}$ is probably due to the fact that the series of fractures parallel to z help in connecting other ones oriented parallel to y.

Two general remarks can be made to conclude this Section. They will also make a useful transition to the next Section.

The macroscopic permeability of fracture networks depends on many geometrical parameters, including the size, shape and orientation distributions of the fractures, and in the present case on the existence of regularly spaced parallel series. A full description should account for all of these parameters. However, it appears that a great part of their influences can be summarized by the dimensionless density $\rho'$, as was already shown by Koudina *et al.* (1998) for isotropic networks, as indicated by the good agreement observed in Fig. 11 between the present calculations and the result (25) of Koudina *et al.* (1998).

Even though Snow's model (29) fails to account accurately for the effect of size polydispersity and of the presence of series of parallel fractures, it still provides a reasonable estimate of the network permeability, considering its low computational cost, including the influence of non-uniform fracture permeability. Of particular interest is the fact that, aside from the fracture permeability, it only requires a few geometrical parameters which are readily available
from line surveys, without need for any assumption regarding the size and shape of the fractures.

# 4 VALIDATION AND EXTENSIONS

This Section is devoted to the validation and the extension of the previous results.

## 4.1 Validation

There are not many examples in the literature where the fracture network is perfectly known and where the previous technique can be validated. However, recently we have studied a sample of fractured granite (Gonzalez Garcia et al. 2000).

This study was conducted on a block of dark grey Hercynian granite from La Peyratte, Deux-Sèvres, France. It is fine-grained (1–2 mm long crystals) and is crosscut by numerous fractures surrounded by discoloured alteration haloes. The primary acquisition was undertaken by Ledésert et al. (1993). The granite block (about 52 × 35 × 36 cm³) was sawed into nine parallel plates, 4 cm in thickness. Trace maps were drawn from the nine sections, by visual examination of the alteration zones, due to the circulation of hydrothermal fluids (Fig. 12). Three examples of these traces are given in Fig. 12. The fracture pattern appeared to be composed of two main families A and B, at about ±30° inclination angle from the vertical axis. These two sets are associated with one horizontal fracture (fracture number 16, in Fig. 12).

This block and its properties were extensively studied by Gonzalez Garcia et al. (2000). For our purposes here, this paper can be summarized as follows. The block of volume $V = 0.066 \text{ m}^{-3}$ contains $N_f = 90$ fractures which are present in more than one plane. This network possesses $N_l = 214$ fracture intersections. Two basic quantities can be derived from these numbers, namely the volumetric fracture density $\rho$ and the dimensionless fracture density $\rho'$:

$$ \rho = \frac{N_f}{V} = 1364 \text{ m}^{-3}, \quad \rho' = \frac{2N_f}{N_f} = 4.76. \qquad (42) $$

Moreover, the total fracture area is equal to $S = 2.05 \text{ m}^2$; an estimation of the lateral extension of the fractures can be made by assuming that they are equal discs of radius $R$. Then, $S = N_f \pi R^2$, and $R$ is found to be equal to 0.085 m.

The transport properties of the block were also analysed. Using the methodology explained in Section 2.6, the Darcy equation was solved in the real network. The following dimensionless results were obtained along the three axes:

$$ K_{nu,xx} = 0.66, \quad K_{nu,yy} = 1.89, \quad K_{nu,zz} = 0.26, \qquad (43) $$

where the permeabilities are made dimensionless by the factor $\sigma/R$, where $R$ is the previous value of the radius when the fractures are assumed to be circular.

Let us now apply the line analysis detailed in Section 2.1 to this granite sample. More precisely, as seen in Fig. 12, we have selected three of the nine planes and measured the number of intersections either along the horizontal or the vertical axis: in the former case, one has a total of $n_h = 44$ intersections for a total length $L_h = 3 \times 0.52 \text{ m} = 1.56 \text{ m}$; in the latter case, one has $n_v = 18$ intersections for a total length $L_v = 3 \times 0.35 \text{ m} = 1.05 \text{ m}$.

For the sake of simplicity, we shall now assume that the fracture network is made of an isotropic distribution of monodisperse discs of radius $R$ with a volumetric density $\rho$. (4) can be generalized in a straightforward way. The total number $n_t$ of intersections between

<figure><img src="image_16.png" /><figcaption>Figure 12. Three of the nine successive trace maps. The traces of two fractures (1 and 16) are indicated in each section. (Reprinted with permission from Ledésert et al. 1993).</figcaption></figure>

a line of length $L$ and this network is given by

$$ n_t = \frac{\pi}{2} \rho L R^2. \qquad (44) $$

Since the excluded volume of discs is equal to $\pi^2 R^3$, the dimensionless density $\rho'$ can be derived as

$$ \rho' = 2\pi n_t \frac{R}{L} = 2\pi(n_h + n_v) \frac{R}{L_h + L_v}. \qquad (45) $$

It is straightforward to make a numerical application of this formula:

$$ \rho' = 12.7. \qquad (46) $$

This quantity should of course be compared to (42). It can also be compared, and with a better agreement to the prediction $\rho' = \pi RS/V = 8.3$, based on the global fracture area, but the latter requires a knowledge that can only be obtained from 2-D or 3-D data.

Let us evaluate the dimensionless permeability of this equivalent network. We can either use the percolation, like formula (25), or the isotropic Snow formula (28) expressed for discs to obtain the two
following evaluations denoted by $K'_p$ and $K'_s$:

$$K'_{p} = 1.9, K_{s} = 2.7, \qquad (47)$$

values that have to be compared to the average one-and-a-half of the three values given by (43).

Let us now summarize and critically evaluate our findings. It is seen that the dimensionless density is overestimated by a factor slightly smaller than 3. Despite this difference, the permeability estimations are excellent since permeability is at most overestimated by a factor 2.

The difference between (42) and (46) is at least partly due to truncation effects. As seen in Fig.12, many fracture traces are truncated by the sample boundaries, and quite a few of them probably intersect out of these boundaries. These intersections are not accounted for in $N_l$, and therefore, $\rho'$ in (42) is underestimated. Similarly, some fracture intersections within the sample have not been detected because of the relatively large spacing of the section planes. These effects have a much smaller influence on the intersection count along a profile which is basically a measurement of the odometric surface area and thus on the prediction of Snow's formula (28).

During the course of these comparisons, one should have in mind the considerable experimental and numerical labor necessary in Gonzalez Garcia et al. (2000) and the straightforward character of the application of the line analysis which yields results in no time within a factor 3 for the dimensionless density and within a factor 2 for permeability.

The comparison could be made more precise by using formulae such as (29).

## 4.2 Extension to variable surface permeabilities

Let us come back to the problem of determining the surface permeability $\sigma$ of the fractures. Note that the flow behaviour due to the complexity of a fracture's microscopic geometry, is not the focus of the present study. For practical purposes, we can summarize previous works by Mourzenko et al. (1995, 2001). If the fractures are viewed as homogeneous above some length, three major quantities are necessary to characterize them. The first one is the fracture aperture $b_o$ as already used in (20). In order to be more precise, one needs to know more about the statistical organization of the fractures. Generally speaking, a fracture can be considered as two rough surfaces put one on top of one another: the heights $Z(r)$ of these surfaces are distributed according to a Gaussian law with zero mean and variance $\sigma_h^2$. Another statistical ingredient is the correlated character, or not, of the two surfaces. Finally, one needs to characterize the statistical organization within each fracture: this is usually done by the covariance $C_Z(\mathbf{u})$ of the heights within the surfaces which is defined as

$$C_Z(\mathbf{u}) = \langle Z(r)Z(r+\mathbf{u}) \rangle, \qquad (48)$$

where the brackets denote the statistical average over the surface. Very often, this covariance can be characterized by a Gaussian variation:

$$C_Z(\mathbf{u}) = \sigma_h^2 \exp \left[ -\frac{\mathbf{u}^2}{l^2} \right], \qquad (49)$$

where $l$ is a characteristic length.

Master curves were given by Mourzenko et al. (1995) in order to estimate more precisely the fracture permeability as functions of $b_o$, $\sigma_h$ and $l$. Of course, this more precise estimation requires more knowledge on the structure of the fractures than the simple estimate (20).

It may also happen that the fracture characteristics depend on its orientation. This can occur when the apertures depend for instance on some regional orientation of the external stresses. Let us give a few examples of how the previous approaches can be generalized to such a case.

First, it should be noticed that the inclusion of a variable fracture permeability does not present any difficulty in the numerical approach presented by Koudina et al. (1998). Each fracture is meshed by triangles whose size is smaller than a typical size $\delta$, which is chosen by the user; typically, this size is smaller than $R/5$ if $R$ is an order of magnitude of the lateral extension of the fractures. A different permeability can be assigned to each triangle.

Secondly, it should be noticed that the Snow formula (27) is valid when the surface permeability depends on the orientation $\mathbf{n}$. It is also possible to include the case where fractures with the same orientation have different permeabilities.

Thirdly, a whole series of calculations is conducted when the fracture permeability is set equal to unity for the events in the E-W direction and to 1/2 for the others. This specific study was motivated by the geological setting as detailed in Section 3.1.

In this situation, (cases V in Tables 3b and c), $\overline{K}_{Sn}$ decreases and the anisotropy is modified, with an increase of the permeability along the y-axis, which corresponds to the E-W direction.

Moreover, consider the pairs of cases in lines (3, 4) and (7, 8) in Table 3(b). They differ only by the fracture permeabilities, which are either constant and set equal to unity (case U), or set equal to unity for the events in the E-W direction and to 1/2 for the others (case V). The permeability is of course smaller in the latter case, since $\sigma$ is equal or smaller than in case U, whereas the network geometry is identical. Moreover, the permeability anisotropy in the y direction, i.e. the direction of the fractures with a large permeability, is further increased. It can be noted that the model (29) correctly predicts the decrease in mean permeability and the change in anisotropy due to the variations of $\sigma$, although it failed to account for the effect of geometrical polydispersivity.

Finally, the four last lines in Table 3(c) compare bidispersed samples, with constant or variable fracture permeability, and $n_{max} = 0$ or 5. They give rise to the same remarks as the other previous tests. Since the connectivity is intermediate ($\rho' = 7.6$ for $n_{max} = 0$), $\overline{K}$ increases by a larger factor (about 1.8) than the number of fractures $N_f$ (4/3) or the density $\rho'(3/2)$, whereas $\overline{K}_{Sn}$ increases by about 3/2. Again, the $x$ direction becomes the least favourable to the flow. However, the permeability anisotropy is now shifted towards the z-axis, instead of the y-axis. The prediction (29) is in fair agreement with this trend. For instance, the shapes of the tensors $\overline{K}$ and $\overline{K}_{Sn}$ in the last line of Table 3(c), which corresponds to bidisperse networks with non-uniform fracture permeabilities and series of parallel fractures for some of the events, are in remarkable agreement.

In order to conclude this Subsection, it can be said that any extra information can be used without any difficulty to estimate more precisely the permeability of the fracture networks.

## 4.3 Extension to polydisperse networks

This is obviously an important issue since it has been seen in many circumstances that the observed networks have sizes which are distributed according to power laws. This important property was first noticed by Segall & Pollard (1983) and it has been confirmed in many subsequent works. The probability density $h(R)$ of the fracture radius $R$ is given by

$$h(R) = \alpha R^{-\alpha}, \qquad (50)$$

where the exponent $\alpha$ is usually ranging between 1 and 3.
In principle, it is contradictory to utilize a complex fracture size distribution when measurements are only performed along a line. This is the basic reason why we did not want to extend the numerical applications too much in this direction; hence, we restricted ourselves in Section 3 to bidisperse networks which represent only a small complication with respect to the monodisperse case; moreover, this bidisperse character is related to the general properties of the fractures in this region which is characterized by discontinuities running east-west.

However, it may be that one has extra information on a given fracture field which has only been measured along a line: this information, which can come from an outcrop where trace distributions have been measured (Berkowitz & Adler 1998), may indicate a strong polydispersity that may or may not be symbolized by the power law (50). This information can be utilized to make the previous predictions more realistic.

First, it is not difficult to generalize the analysis made in Section 2.2. (4) is valid for fractures of a given direction and a given radius. If one considers fractures with a given unit normal $n$, one has to simply add up the contributions of the various sizes weighted by the probability density (50).

Secondly, it is not difficult to generalize the numerical calculations relative to the percolation threshold and the network permeability for this situation. It should be mentioned here that such calculations are currently made in a systematic way and that general results and master curves will be soon published.

For the time being, one can provide the reader with a simple extension of the Snow formula (28). If the fracture sizes follow the power law (50), the total volumetric area $S$ is given by

$$S = \int_{R_m}^{R_M} \rho A_p(R)h(R) dR = \rho \langle A_p \rangle, \quad (51)$$

where $R_m$ and $R_M$ are the minimal and maximal radii of the discs, respectively. $A_p$ is equal to $\pi R^2$.

This expression can be introduced into (27) and it yields

$$K_{Sn}^{iso} = \frac{2}{3} \rho \langle A_p \rangle. \quad (52)$$

It is also possible to account in (51, 52) for a correlation between the size $R$ and the permeability $\sigma$ of the individual fractures:

$$K_{Sn}^{iso} = \frac{2}{3} \rho \langle A_p \rangle, \langle A_p \rangle = \int_{R_m}^{R_M} \sigma(R) A_p(R) h(R) dR. \quad (53)$$

# 5 CONCLUDING REMARKS

Macroscopic properties such as the volumetric densities and the percolation properties of a fracture network, can be estimated from a line survey provided that the orientation of the fractures are recorded. One of the major merits of these estimations is to provide analytical expressions where the measured and assumed quantities are clearly distinguished. For instance, the percolation character (or not) of the network can be estimated if the lateral extensions of the fractures are known: such quantities cannot be derived from the line data themselves, but can be obtained from other sources such as measurements on an outcrop.

Another important quantity which can be estimated is the macroscopic permeability; of course, this necessitates the creation of an assumption (or possibly measurements or estimates) of the individual fracture permeability and of the lateral extensions of the frac-

tures. Whatever the formula used (see 25 and 29), it can be entirely expressed in terms of measured quantities.

This approach has been validated on a granite block whose structure is known. Its connectivity and permeability properties are satisfactorily estimated by the proposed approach. One can also appreciate with this example, the drastic reduction in the amount of work required to fully calculate these properties, compared to the one necessary for the present technique.

Several extensions to variable fracture permeabilities and polydisperse networks are also proposed.

This work will be continued in several ways. The properties of these anisotropic polydisperse networks composed of a series of parallel fractures require a complete study. Moreover, the initial motivation of this study, which was the evolution of these karstic fields under the influence of dissolution, will be completed.

# ACKNOWLEDGMENTS

This work has been partly supported by the INSU programme PNRH and by the European contract Saltrans EVKI-CT-2000-00062. Most calculations were done at CINES (Montpellier). The authors also wish to thank E. Debros for kindly providing unpublished data relative to the geology of the Baget watershed (cf. Fig. 2).

# REFERENCES

* Adler, P.M., 1992. Porous Media: Geometry and Transports, Butterworth/Heinemann, Stoneham, MA.

* Adler, P.M. & Thovert, J.-F., 1999. Fractures and Fracture Networks, Kluwer, Dordrecht.

* Balberg, I., Anderson, C.H., Alexander, S. & Wagner, N., 1984. Excluded volume and its relation to the onset of percolation, Phys. Rev., B30, 3933-3943.

* Berkowitz, B. & Adler, P.M., 1998. Stereological analysis of fracture network structure in geological formations, J. geophys. Res., B103, 15 339-15 360.

* Charlaix, E., 1986. Percolation threshold of a random array of discs: a numerical simulation, J. Phys., A19, L533-L536.

* Choukroune, P., 1992. Tectonic evolution of the Pyrenees, Ann. Rev. Earth planet. Sci., 20, 143-158.

* Debros, E.J., 1987. Modèle de bassin triangulaire à l'intersection de décrochements divergents pour le fossé alco-cénomanien de la Ballongue (zone nord pyrénéenne, France), Bull. Soc. Geol. Fr., 5, 887-898.

* Fisher, M.E., 1971. The theory of critical point singularities, in, Critical phenomena, Proc. 51st Fermi School, Varenna, Italy, ed. Green, M.S., pp. 1-99, Academic Press, New York.

* Gonzalez Garcia, R., Huseby, O., Thovert, J.-F., Ledésert, B. & Adler, P.M., 2000. Three-dimensional characterization of fractured granite and transport properties, J. geophys. Res., 105B, 21 387-21 401.

* Goula, X., Olivera, C., Fleta, J., Grellet, B., Lindo, R., Rivera, L., Cisternas, A. & Carbon, D., 1999. Present and recent stress regime in the Pyrenees, Tectonophysics, 308, 487-502.

* Huseby, O., Thovert, J.-F. & Adler, P.M., 1997. Geometry and topology of fracture systems, J. Phys., A30, 1415-1444.

* Johnson, J.A. & Clarence, A.H., 1989. The structural and sedimentary evolution of the Cretaceous North Pyrenean Basin, southern France, Geol. Soc. Am. Bull., 101, 231-247.

* Kingman, J.F.C., 1993. Poisson Processes, Oxford Science Publications, Oxford.

* Koudina, N., Gonzalez Garcia, R., Thovert, J.-F. & Adler, P.M., 1998. Permeability of three-dimensional fracture networks, Phys. Rev., E57, 4466-4479.

* Ledésert, B., Dubois, J., Velde, B., Meunier, A., Genter, A. & Badri, A., 1993. Geometrical and fractal analysis of a three-dimensional hydrothermal vein network in a fractured granite, J. Volc. Geotherm. Res., 56, 267-280.
Long, J.C.S. & Billaux, D.M., 1987. From field data to fracture network modeling: an exemple incorporating spatial structure, *Water Resources. Res.*, **23**, 1201–1216.

Mangin, A., 1974. Contribution à l'étude hydrodynamique des aquifères karstiques, *Ann. Splol.*, **29**, 495-601.

Mourzenko, V.V., Thovert, J.-F. & Adler, P.M., 1995. Permeability of a single fracture; validity of the Reynolds equation, *J. Phys. II France*, **5**, 465-482.

Mourzenko, V.V., Thovert, J.-F. & Adler, P.M., 2001. Permeability of self-affine fractures, *Transp. Porous Media*, **45**, 89-103.

Pistre, S., Lopez-Chicano, M., Pulido-Bosch, A. & Drogue, C., 1999. The role of western Mediterranea tectonic evolution in the geometry of a karstic domain in the Betic Cordilleras (Sierra Gorda, Spain): importance of a tardy extensional regime, *Geodin. Acta*, **12**, 11-24.

Reinecker, J., Heidbach, O. & Mueller, B., 2003. *The 2003 release of the World Stress Map*, available online at www.world-stress-map.org

Reynolds, P.J., Stanley, H.E. & Klein, W., 1980. Large-cell Monte Carlo renormalization group for percolation, *Phys. Rev.*, **B21**, 1223-1245.

Segall, P. & Pollard, D.D., 1983. Joint formation in granitic rock of the Sierra Nevada, *Geol. Soc. Am. Bull.*, **94**, 563-575.

Snow, D.T., 1969. Anisotropic permeability of fractured media, *Water Re-sour. Res.*, **5**, 1273-1289.

Souriau, A. & Granet, M., 1995. A tomographic study of the lithosphere beneath the Pyrenees from local and teleseismic data, *J. geophys. Res.*, **100**, 18 117-18 134.

Souriau, A., Sylvander, M., Rigo, A., Fels, J.F., Douchain, J.M. & Ponsolles, C., 2001. Sismotectonique des Pyrnes: principales contraintes sismologiques, *Bull. Soc. Geol. Fr.*, **172**, 25-39.

Stauffer, D. & Aharony, A., 1994. *Introduction to Percolation Theory*, 2nd edn, Taylor and Francis, Bristol.