## Abstract

Field observations from several outcrops in the Eastern High Atlas Mountains, near Amellago (Morocco), are used to determine fracture-network model parameters, such as the aspect ratio of the fractures represented as rectangles whose longer side is horizontal, the volumetric area of fracture surfaces, the fracture mean size and the fracture density. The fracture orientations can be roughly approximated by Fisher distributions, where the parameters are determined by outcrop measurements. The permeability of the fracture networks can be calculated by application of the Snow equation for infinite fractures or by numerical resolution of the flow equation for fracture networks generated with the parameters deduced from the outcrop measurements. These two permeability estimations are shown to be in good agreement, which suggests that theoretical or semi-empirical solutions may provide reasonable approximations of fracture-network permeability in some carbonate reservoirs when conditioned to appropriate outcrop and subsurface data.

## Introduction

Fractures influence the hydraulic and transport properties of natural geological formations. The purpose of the present work is to predict the properties of a real reservoir comprising fractured carbonates; in particular, its percolation status and its permeability.

These properties can be obtained either by using theoretical predictions for some standard cases, such as the Snow model (Snow 1969) for very well connected fracture networks, or by applying direct numerical simulations. Both approaches are based on statistical geometrical characteristics of the fracture networks. In the first case, existing models are directly applied. In the second case, the hydraulic properties are obtained by solving the governing flow equations, which are written on the local scale and applied to the fracture network represented in an explicit way. Since a real fully described, three-dimensional (3D) sample cannot be provided, reconstructed ones are obtained based on the statistical geometrical characteristics of the real networks. Therefore, these characteristics are essential for the prediction of the hydraulic properties, whatever the method.

The statistical geometrical characteristics are extracted from outcrop data. Commonly, information describing such characteristics is limited to the intersection counts along boreholes or the trace maps on outcrops. Several stereological methods were used to determine the required parameters from the available field measurements (Berkowitz & Adler 1998; Sisavath *et al.* 2004; Thovert & Adler 2004; Gupta & Adler 2006; Mourzenko *et al.* 2011*a*).

The present paper is organized as follows. The geological background and the theoretical framework are given first. Then, the fracture-network model parameters and a few hypotheses are described. The observations from four vertical cliff-face outcrops and one pavement outcrop in the Eastern High Atlas Mountains, near Amellago (Morocco), and the corresponding trace maps are analysed. The treatment and subsequent analysis of the available field data enable the determination or calibration of the model parameters. Then, the permeability tensor is derived from Snow’s formula, which is valid for infinite fractures, or equivalently for high fracture density (Mourzenko *et al.* 2011*a*), and compared to the results of direct calculations in numerical samples of fracture networks generated according to the model parameters.

### Geological Background

The primary focus of this paper is to demonstrate a method for rapid estimation of fracture permeability. We include a brief overview of the geological framework here as context for the natural-fracture data used as inputs for the permeability modelling. The overview draws, in part, from publications related to studies of the same outcrop that formed part of a large multidisciplinary research effort (Agar *et al.* 2010; Amour *et al.* 2012, 2013; Christ *et al.* 2012; Shekhar *et al.* 2014). The study location is situated in the eastern portion of the High Atlas Mountains where Jurassic sediments have been exposed by uplift related to inversion along the ENE–WSW-trending margin of a Triassic–Jurassic rift system and subsequent erosion (Jacobshagen *et al.* 1988; Warme 1988; Beauchamp *et al.* 1996) (Fig. 1). Triassic rifting in the Atlas region was localized along a WSW–ENE-trending fault system and connected to the Western Tethys. The rift basin separated the Saharan craton (Gondwana) from two microplates, the Moroccan and Oran Mesetas (Warme 1988). Several kilometres of continental sediments accumulated within this rift through the Triassic until the Early Jurassic (Warme 1988). Jurassic transgression was marked by a shift to a marine depositional environment in the central and eastern High Atlas rift, generating large carbonate ramp systems (Christ *et al.* 2012).

The carbonate sediments that host the fractures discussed in this paper originated within a Jurassic system of two carbonate platforms (Lias and Dogger in age) separated by a 100 m-thick Toarcian hemi-pelagic marl succession (Du Dresnay 1979; Warme 1988). Amour *et al.* (2012) suggested that the drowning of the older carbonate platform might have resulted from a renewed phase of extension and/or a sea-level rise with an associated Toarcian oceanic anoxic event. The shallow-water deposits of the upper carbonate platform, observed in this study, prograded to the NE (Stanley 1981). A transtensional tectonic regime that formed numerous rhomb-shaped basins has been interpreted from the Toarcian to the Bathonian (Du Dresnay 1979; Crevello 1990). It has been postulated that the pull-apart structures controlled local depositional relief, with shallow-water carbonates accumulating on local highs associated with structurally tilted blocks (Stanley 1981). Further extension during the Aalenian–Early Bajocian drove a broad change from a rimmed platform to a carbonate ramp to the east of the study area (east of Rich: Fig. 1). The region surrounding the outcrop discussed in this study remained stable, with strata recording two progradations of an oolitic carbonate ramp towards the basin, now referred to as the Amellago Formation (Late Toarcian–Early Bajocian age) and the Assoul Formation (Middle–Late Bajocian age) (Hadri 1993; Durlet *et al.* 2001; Pierre *et al.* 2010). The Agoudin marls of Early Bajocian age separate the deposits of the two ramp systems (Fig. 2).

With the Cenozoic convergence of Africa and Europe, the Jurassic strata were translated and uplifted along numerous reverse and thrust faults, converging mainly to the south. Although not mapped in detail, the hemi-pelagic marls appear to have provided numerous weak zones that accommodated contraction during subsequent inversion. The contractional deformation of the eastern High Atlas is more intense towards the frontal thrust zones. A transect from Goulmina (south) to Amellago (north) (Fig. 1) passes through numerous zones of intense thrust-related deformation (100–500 m wide) and penetratively deformed strata, with exposures of tight to isoclinal folds. Closer to Amellago and further north, individual thrust zones become more widely spaced, with deformation mainly localized to fault damage zones (*c*. 0.5–20 m wide). With these changes, the fold morphologies also change to dominantly broad flexures (km-scale wavelengths) and, occasional, locally tight flexures immediately adjacent to reverse fault planes (mainly truncated folds in fault hanging walls).

The outcrop focus of this study is located in the Amellago Canyon, which is located approximately 65 km to the SW of the town of Rich and 8 km NW of the village of Amellago. The outcrop is referred to here as the ‘Island’ based on its isolated relief, formed by fluvial erosion of the Gheris river (Fig. 1), that provides approximately 90% exposure of the strata. The Amellago Canyon provides tremendous exposures of carbonate strata in cliffs reaching well over 150 m in places. The elliptical Island (*c*. 320×190 m) offers exposures of strata from the Assoul Formation from all directions. The Assoul Formation is approximately 300 m thick and comprises shallow-water deposits that prograded towards the NE sub-basin (Poisson *et al.* 1998, Pierre *et al.* 2010; Amour *et al.* 2012; Christ *et al.* 2012). It comprises alternating deposits of two types of carbonate ramp, an ooid-free muddy ramp and an oolitic ramp, which formed during the early transgressive and late transgressive–highstand phases of fourth-order depositional sequences, respectively (Pierre *et al.* 2010; Christ *et al.* 2012). Both of these ramp types contain bioclastic and peloidal wackestone–packstone to floatstone from the inner to proximal middle ramp, with local occurrence of coral-dominated bioconstructions.

The strata exposed on the Island dip gently, overall, to the north (*c*. 12^{o}). They have been gently flexed by contractional deformation, interpreted as being related to the Atlas inversion (Fig. 3). The Island stratigraphy (Amour *et al.* 2012; Christ *et al.* 2012), combined with striations and polyphase slickenside development on east–west-trending fault surfaces, indicate the presence of low-offset (*c*. 1–10 m) fault zones within the Island and the surrounding area. The two most distinct fault zones cut across the southern end and centre of the Island, striking east–west (Fig. 3). The steeply dipping southern fault zone locally reaches over 10 m in width and contains brecciated carbonate with polyphase carbonate cements and minor pyrite mineralization. On the Island, there is an apparent normal offset of approximately 10 m down to the south (based on the stratigraphic interpretation of Amour *et al.* 2013). Striations and slickensides within the fault zone support local oblique movements but we cannot be sure that these formed synchronously with the movements that generated the apparent normal offset. In contrast, the fault zone in the centre of the Island exposes minimal fracture damage (a few centimetres at most) around a discrete, steeply south-dipping fault plane with patchy carbonate cements preserved. The fault has a normal displacement component of approximately 1 m down to the south.

In addition to the low-offset faults, the Island is dissected by penetrative fracture populations and abundant bed-parallel stylolites with occasional moderate to steeply dipping stylolites. Some fractures are bed-confined or terminate at stylolites within beds, while others dissect multiple beds, extending vertically for over 10 m in places. Field observations support the existence of four sets of fractures striking east–west, ESE-WNW, ENE-WSW and north–south (weakly developed). Most fractures have apertures of the order of 1 mm or less. The dominant fracture fill is calcite with minor clays, but many of the exposed fractures contain either no cement or sparse patches of cement. In many cases, we were unable to determine whether cements have been eroded or no cements ever formed. Clear, cross-cutting relationships among different fracture sets are rare, preventing any rigorous interpretation of the relative timing and evolution of the fracture system. We observed a few fractures on the Island that were overprinted by north-dipping stylolites, suggesting that at least some of the fractures formed during or prior to contraction (similar relationships have also been observed in other sections of the Amellago Canyon: Agar *et al.* 2010). Based on the overall structural framework, it is possible that the north–south-trending fractures formed in response to compressional stresses driving the north–south Atlas inversion, while the east-, ESE- and ENE-trending fractures may be related to weak flexure during inversion and late oblique-normal and oblique-reverse faulting. The absence of cements in many fractures may indicate relatively late origins that could be related, at least in part, to stress release during exhumation. For the purposes of this study, we are concerned with the final cumulative population of fractures. The uncertainty in the relative timing of fracture formation is not significant for the method we aim to demonstrate.

The data used in this study come from four ‘windows’ studied on the Island (Fig. 4). Each window is approximately 30 m long and 5 m high. The selected height of the windows was controlled, to a large extent, by the weathering profile of the outcropping strata. Fractures were measured in resistant strata that form small cliffs circumscribing the Island, separated by intervals of weathered marls. Most of the measured fractures terminate within or at the top or bottom of each cliff. In a few cases, it was not possible to confirm the nature of the fracture terminations owing to weathering in strongly fissile marls above or below the cliff. In the case of windows 4 and 6, where taller cliffs were present, the longer vertical extent of a few widely spaced fractures was documented. Measured fracture orientations were tied to line drawings of fractures constructed for each window. The line drawings show fracture traces on photomontages that have been corrected for perspective and distortion associated with the camera lens and image acquisition from locations at slightly topographically lower levels than the base of the window. The following descriptions outline the characteristics of each window.

Window 4 is located at the topographically lowest position on the eastern side of the Island, close to the lowest extent of rock exposure. It comprises skeletal wackestone–packstone interlayered with bioclastic mudstone and marl intervals, and shoals of cross-bedded ooid grainstones. There are local build-ups or ‘mud mounds’ with oysters and coral fragments, which are interpreted to have formed in a middle ramp setting.

Window 5 is located within a small-scale shallowing-upwards succession within a medium-scale deepening-upwards succession on the western side of the Island, approximately 30 m below the highest point on the outcrop and directly below a hardground (Amour *et al.* 2012). It comprises mainly bioclastic peloidal wackestone overlying marly oyster beds in a proximal middle-ramp setting. Some of the fractures in this window have very small (1–10 cm) apparent normal offsets and appear to accommodate early, very minor extension, with displacements dying out into the marl beds at the base of the exposure.

Window 6 is located above a hardground (HG4 of Amour *et al.* 2012) on the eastern side of the Island, approximately 20 m below the top of the outcrop. It comprises oncoidal floatstone, oncoidal mudstone and wackestone–packstones with oncoids, and is interpreted to have formed in an inner-ramp setting, representing a continuation of the small-scale shallowing-upwards trend from Window 5.

Window 7 is situated just below a horizon interpreted to be a maximum flooding surface at the top of a small-scale deepening-upwards succession within a medium-scale shallowing-upwards succession. It comprises mainly bioclastic wackestone and bioclastic packstone (gastropods, corals, oyster fragments) with local marls, and a thin (<0.25 m) interval of bioclastic floatstone at the top of the window. These sediments are interpreted to have formed in a proximal middle-ramp setting. The south end of this window terminates close to the low-offset normal fault described above. There is no discernible increase in mesoscopic fracture intensity towards the fault.

The analysis presented here also uses fracture trace-length data from a pavement exposed at the southern end of the Island. The pavement is located to the south of ‘F2’ (Fig. 3a) and formed on top of strata that represent the lateral continuation of Window 7. The eroded planar surface exposes the horizontal traces of fractures dissecting the Island area and is approximately 35 m long (north–south), with a maximum width of approximately 25 m (east–west). In several cases, the measured fracture trace lengths do not capture the maximum lengths due to truncation by the eroded pavement edge or burial beneath rubble that has accumulated in places on the pavement.

## Theoretical Framework

Most symbols are defined in Table 1.

### Geometrical characteristics of fracture networks

Let us consider a 3D network of plane, convex and, possibly, polydisperse fractures, which is characterized by the mean fracture area,
, the mean fracture perimeter,
, and the density, *ρ*, defined as the number of fractures per unit volume. The brackets denote the statistical average. The fractures are assumed to be convex, with an identical shape characterized by the shape factor, *η*_{P}:
(1)

where *R* is the radius of the disk circumscribed to the fractures.

The fractures are isotropically orientated within their plane. The normals to the fracture ** n** are usually expressed in the standard polar coordinates

*θ*and

*ψ*illustrated in Figure 5. The probability density function (pdf) of

**is denoted by**

*n**f*(

*θ, ψ*). A classical form is the Fisher distribution characterized by the pole

**. When**

*p***is given by its polar coordinates (**

*p**θ*

_{p},

*ψ*

_{p}),

*f*(

*θ, ψ*) is written as: (2)

where *κ* is a parameter that tends towards zero when the distribution tends to be isotropic; when *κ*→∞, all of the normals are parallel to ** p**.

However, the fracture orientations can be provided by a set of field observations. Such a dataset can be used directly, or modelled by a Fisher distribution, or the superposition of several families obeying Fisher distributions in order to apply generic models whenever possible.

The fracture density *ρ* is the number of fractures per unit volume. It is useful to introduce the dimensionless density, *ρ*′, defined as the average number of intersections per fracture, which can be expressed as:
(3a)

where *γ* is the angle between two fractures. This expression makes use of the concept of excluded volume, as introduced in this field by Balberg *et al.* (1984). If the fractures are isotropically orientated,
and the dimensionless density is given by (Adler *et al.* 2012):
(3b)

For anisotropic networks, *ρ*′, can be expressed as:
(3c)

where Φ is a correction factor, . For a Fisher distribution, Φ is given by: (4)

where *I*_{0} and *I*_{1} are the two first modified Bessel functions.

### Stereological analysis

Some geometrical characteristics of fracture networks can be directly related to measured 1D or 2D parameters. The intersections of the fractures with a plane are called traces or chords.

The average number of intersections per unit length,
, of traces with a line parallel to a unit vector ** p** can be easily measured; it is equal to:
(5a)

When there is no correlation between the fracture size and orientation, equation (5a) can be written as: (5b)

Then, 2D measurements are considered. The average number of traces per unit area *Σ*_{t} is given by Mourzenko *et al.* (2011*a*) as:
(6a)

where α is an angle between the normal *n*_{w} to the observation plane and the normal **n** to the fracture. When there is no correlation between the fracture size and orientation, equation (6a) can be written as:
(6b)

The trace length per unit area *C* is given by:
(7a)

When the fracture orientations and sizes are independent, equation (7a) reads as: (7b)

The average trace length, , is measured as the total trace length divided by the number of traces: (8)

The general relationships (equations 5b, 6 & 7) between the measured quantities and the 3D geometrical characteristics were determined for Fisher distributions (Mourzenko *et al.* 2011*a*).

### Hydraulic properties of fracture networks in generic cases

A fracture network percolates when it extends from − ∞ to + ∞ in at least one direction of space. For monodisperse isotropic networks, the percolation threshold (*ρ*′_{c}) above which networks percolate is found to be nearly independent of the fracture shape and equal to 2.3 (Huseby *et al.* 1997). A correction is necessary only for very elongated fracture shapes, in which case Mourzenko *et al.* (2005) found that:
(9)

Mourzenko *et al.* (2009) showed that *ρ*′_{c} is not significantly modified when the fracture planes are orientated according to a Fisher distribution when *κ* is ranging between 0 and at least 50.

If the network percolates, its permeability is obtained as follows. At a scale intermediate between the typical fracture aperture, *b*_{f}, and its overall extent, the flow inside each fracture is governed by the 2D Darcy law and the mass conservation equation:
(10)

where ** q** is the locally averaged flow rate per unit width,

*μ*is the viscosity, σ is the fracture transmissivity, is the 2D gradient operator in the mean fracture plane and is the pressure gradient in this plane. In the present study, σ is taken to be uniform over each fracture and identical for all the fractures.

When a macroscopic pressure gradient
is applied to the fracture network contained in the volume τ_{0}, the seepage velocity
of the resulting flow is evaluated by:
(11a)

where *τ*_{f} is the interstitial volume of the fractures and *S*_{f} their projection onto the mean planes. The flux is related to the pressure gradient by Darcy’s law:
(11b)

where ** K** is the fracture network permeability.

Usually, ** K** is derived numerically, but it can be derived analytically for a network of infinite fractures with an average surface area

*S*of fractures per unit volume. When the fracture transmissivity, σ , is uniform over each fracture and constant for all the fractures, the permeability

**is given by Snow’s formula (Snow 1969) or Oda’s formula (Oda 1985), which was originally derived by Romm (1966): (12a)**

*K*where ** I** is the unit tensor and

**a dyadic product. For isotropic orientations, one obtains: (12b)**

*nn*For anisotropic networks, *K*_{S} can be written as
(13a)

where **Ψ** incorporates the effect of the geometrical anisotropy:
(13b)

For a Fisher distribution with a pole parallel to the *z*-axis, the tensor, **Ψ**, can be written as:
(13c)

where Ψ_{⊥} and Ψ_{||} were found by Mourzenko *et al.* (2011*a*) as:
(13d)

Of course, **Ψ** can also be directly calculated via equation (13b) based on a dataset of orientations measured at outcrop.

Snow’s result (equation 12a) is an upper bound for the permeability of networks of finite fractures. It is a good first approximation when the networks are very well connected; that is, when *ρ*′ is very large. The average surface of fractures per unit volume is given by:
(14)

## Parameters of the Fracture Network and Assumptions

In order to describe fully a fracture network, the following information is required. One should know the fracture positions, orientations and shapes, as well as the size distribution and the density of the fractures.

In order to simplify the model, a few hypotheses about these parameters are used, which are supported by intuitive arguments or by direct observations.

The *fracture positions* are supposed to be random and there is no correlation. This hypothesis is supported by the variograms of the spacings between successive intersections along scanlines, shown later, which do not reveal any spatial correlations.

The *fracture orientations* are directly taken from field measurements. There is no simplifying hypothesis, except for the assumption that the sampling is unbiased in terms of probabilities in numbers. However, a generic modelling by Fisher families is possible; it is not necessary for the reconstruction of fracture networks, but it may be convenient for the analysis, since some theoretical results are available in this case.

The fractures are *rectangular with a constant aspect ratio* and with their long side horizontal or parallel to bedding. The supporting arguments are the following. On the trace maps (e.g. Fig. 4), one can see that the traces often terminate at bedding planes or stylolites, such that they are bounded by the strata. The fractures are elongated; therefore, the exact shape of their tips has relatively little importance, and assigning a rectangular shape to the fractures is the simplest choice. Then, intuitively, fractures with the longest horizontal extent are expected to also have the longest vertical extent. In the absence of data, again, the simplest choice of a constant aspect ratio is made.

There is no correlation between the fracture size and the orientation. The size distribution is described by the number-pdf *n*(*a*), where *a* is the shorter fracture side. It will be demonstrated that the correlation coefficient for the trace lengths and orientations is very small.

The density of the fracture network, *ρ*, is defined as the number of fractures per unit volume. *ρ* is supposed to be uniform over the investigated domain.

For illustrative purposes, the *fracture aperture* is assumed to be uniform in the prediction presented below. This assumption is not as restrictive as it may seem at first sight; Hamzehpour *et al.* (2009) showed that a mean field approximation is in good agreement with detailed calculations performed for heterogeneous fractures.

All of these hypotheses can be criticized since, at best, they represent a simplified reality. The problem is normal for modelling real geological domains. On the one hand, an appropriate numerical model can describe any fracture network (even with non-planar fractures) with different physical properties at each node of the mesh (Adler *et al.* 2012). On the other hand, even for well-documented outcrop datasets such as that from the Amellago Canyon outcrops, geological data are scarce and one has to make assumptions in order to build statistical network models. All of the assumptions made herein cannot be considered as exact, but they provide a representation of the real fracture field, which is broadly consistent with geological interpretations.

## Measurements and Treatments of the Field Data

The observations from the Amellago Canyon outcrops are analysed in this section. For the four vertical outcrop windows described above (Figs 3 & 4), fracture orientations and trace maps drawn from photographs, represented by coordinates of trace end points, were recorded.

Independently, field measurements of fractures on one horizontal pavement (Fig. 3) were also recorded. They comprise a list of fracture trace lengths, the corresponding fracture azimuths and the nature of the trace terminations (inside or outside of the observation domain).

### Three-dimensional fracture orientations

The field measurements provide us with the fracture orientations given by the dip, *θ*_{f}, the dip azimuth, *ψ*_{f}, and the azimuth of each vertical outcrop surface, *ψ*_{0}. The normal vectors, ** n**, to the fractures are expressed in polar coordinates (

*θ*

_{n},

*ψ*

_{n}). The peaks in the orientation distributions are more easily visualized by including both vectors +

**and −**

*n***for each observation. Therefore, the fracture orientations are defined in the whole domain [0°<**

*n**θ*≤180°, 0°<

*ψ*≤ 360°] and the data in different outcrop windows can be compared. In addition, each discrete observation (

*θ*

_{n},

*ψ*

_{n}) is replaced by a distributed density in the form of a Fisher distribution around (

*θ*

_{n},

*ψ*

_{n}): (15)

where *β* is the angle between (*θ*_{n}, *ψ*_{n}) and (*θ, ψ*). This accounts for an uncertainty in the field measurements. The parameter *κ* is set to *κ*=87, which corresponds to a 95% confidence range of 15° around the recorded orientation.

The probability of observing a fracture in an outcrop is proportional to sin α where α is the angle between ** n** and the normal

*n*_{w}of the outcrop surface (Fig. 5b). This is the so-called Terzaghi’s correction (Terzaghi 1965). The real fracture orientation distribution in outcrop

*i*is written as: (16a)

This formula is renormalized in the sense that
is equal to 1. *f*_{4} and *f*_{5} are displayed in Figure 6a, b for the outcrop windows 4 and 5.

The same sets seem to exist in all the studied outcrops, although they are more or less visible depending on the orientation of the outcrop surfaces. This suggests that all of these observations result from a unique orientation distribution sampled in different ways by the various outcrop windows and pavement. This distribution *f*_{tot}(*θ, ψ*) can be obtained by combining all of the data:
(16b)

where *N _{i}* is the number of fractures in outcrop

*i*and . The distribution

*f*

_{tot}(

*θ, ψ*) is displayed in Figure 6c.

*f*

_{tot}(

*θ, ψ*) can be approximated by the superposition of two or four Fisher distributions. Only the second superposition will be described in detail, for the sake of brevity.

Consider the probability distribution function *f*_{tot}(*θ, ψ*) in the half domain [0°<θ≤180°, 52°<ψ≤232°], which eliminates the replication of the data. Four families of normal vector orientations are distinguished, which can be represented by Fisher distributions defined by their poles denoted by *p*_{j,4F} and the corresponding parameters denoted by *κ _{j}*,

_{4F}(

*j*=1, 2, 3, 4). For their determination, the domain is divided into four subdomains as shown in Figure 7b.

In each subdomain *j*, the pole *p*_{j,4F} is found as the direction that minimizes the scatter of the data. More precisely, the tensor **Ψ** (equation 13b) is evaluated by averaging the field data, and diagonalization. The pole *p*_{j,4F} corresponds to the eigenvector of **Ψ** associated with the smallest eigenvalue Ψ_{||}. This is equivalent to maximizing the quadratic mean of the cosine of the angle between *p*_{j,4F} and the fracture normal vectors. Then, *κ _{j}*,

_{4F}is determined by solving equation (13d).

In polar coordinates, the poles and parameters of each family denoted by the superscript (*i*) are obtained as:
(17a)

The total four-peak model distribution is given by: (17b)

where are the fractions of the data that belong to each family: (17c)

The distribution *f*_{4F}(*θ, ψ*), obtained by equation (17b), is a good approximation to *f*_{tot}(*θ, ψ*), with a correlation coefficient around 0.96. The shape of *f*_{4F} in Figure 7d is very similar to the original distribution *f*_{tot} (Fig. 6c).

A similar superposition of two Fisher distributions can be carried out: (18a)

(18b)

(18c)

This superposition is detailed in Figure 7a and its result is given in Figure 7c.

In summary, the field observations from all outcrop windows can be represented by a single orientation distribution, which can be approximated by a two- or four-peak model. When Figure 7c and d are compared, it is apparent that the approximation by four Fisher distributions is much better than the approximation by two.

### Two-dimensional measurements

The first set of 2D data is the end-to-end trace maps in the vertical outcrop windows and the second is the list of trace lengths on the outcrop pavement. Hereafter, the subscripts o and p denote the parameters and measurements in the outcrop windows and in the pavement, respectively. The traces of Window 5 are displayed in Figure 8a, b, as an illustration.

The trace maps are studied in order to measure the trace length per unit area, *C*, the number of traces per unit area, *Σ*_{t}, and the mean trace length,
, which are defined as:
(19)

where *c _{i}* is the length of trace

*i*, Ω is the sampling domain and

*N*

_{t}is the number of trace terminations within the domain.

In each outcrop window, the upper and lower bedding plane or stylolite are used as domain boundaries. The traces that cross the boundary are truncated, and those which are out of the domain are not taken into account.

The domain area, Ω_{o}, the number of trace terminations, *N*_{t,o}, within the domain and their total length are measured. Then, the parameters *C*_{o}, *Σ*_{t,o} and
_{Ω,o} for each outcrop are calculated from equation (19).

Although the domain area of the pavement, Ω_{p}, is not available, the mean trace length
_{Ω,p} can be determined from the trace lengths and the data about the trace terminations.

Consider the statistics of the individual trace lengths recorded in the whole map of each outcrop window and in the pavement. The mean trace length,
_{o,p}, and the corresponding standard deviation, σ_{c}, are measured. The ratio
is found to be nearly the same for all the outcrop windows. This constant, denoted by γ, is determined as:
(20)

Then, the cumulative trace-length distribution in the outcrop windows and in the pavement is well approximated by the following fit, which depends only on the mean trace length
_{o,p}:
(21a)

with (21b)

The probability density function is derived from equation (21a) as: (21c)

In Figure 8c–e, the approximation of the trace-length distributions by Φ(*c*) with parameters *c*_{0} and λ determined by equation (21) is compared with the data, and it is found to be in good agreement with them.

It should be noted that the mean trace length,
, was measured in two different ways that yield slightly different results.
_{o} is obtained by directly averaging the length of the traces recorded in the whole window; that is, by dividing the sum of the trace lengths by the number of observed traces. It probably underestimates
because of truncations by the boundaries of the observation region, which are not precisely known.

When operating on a subdomain Ω in the window, these truncations are even more important; however, since the domain boundaries are explicitly defined, they can be accounted and corrected for.
_{Ω} is defined in equation (19) as the ratio *C*/*Σ*_{t}. The trace length density, *C*, is a robust estimator since it does not depend on the domain size, except for statistical fluctuations. Conversely, the trace number density, *Σ*_{t}, is evaluated as half the density of trace terminations. Thus, it takes into account the truncations by the boundary ∂Ω.

In each outcrop,
_{Ω} is slightly larger than
_{o} and represents the expected value of
_{o} if the outcrop is large enough for truncation effects to be negligible.

A scanline analysis of the full trace-map data was performed. The scanlines are set parallel to bedding. The distribution of the spacings, *s _{i}* (Fig. 9), between the successive intersections of scanlines with the traces is in good agreement with the cumulative distribution function:

(22)

where
is the mean spacing. The distribution *F* in equation (22) is the expected result if the intersections are randomly distributed according to a Poisson process, without any spatial correlation.

This is confirmed by the semivariogram γ(*k*) for the spacings *s _{i}*, defined as:
(23)

In Figure 9, γ(*k*) reaches the variance
as soon as *k* ≥ 1 and remains constant within some small fluctuations. This suggests that there is no correlation between successive spacings, *s _{i}*. Therefore, the hypothesis about random fracture positions in the model is supported by the data.

Finally, it was shown that there is no correlation between the fracture size and its orientation. The intercorrelation coefficient between these two quantities was systematically calculated; for Window 5 and for the horizontal pavement, it is equal to −0.2 and −0.08, respectively. For the sake of brevity, the details are not reported here.

## Geometrical Characteristics of the Fracture Networks

The fractures are assumed to be rectangles, with their longer side orientated horizontally. Let *a* and *b*(*a*) be the shorter and longer sides of the fracture, respectively, *f* is the constant aspect ratio, *b*(*a*)=*fa*, and *ρ*(*a*) is the volumetric density. It is assumed that there is no correlation between the fracture length and the orientation. Then, the parameters *Σ*_{t,o} and *C _{o}*, which are defined by equation (19), can be evaluated for each outcrop window as:
(24a)

(24b)

where *ρ*(*a*) is the volumetric density, which is related to the fracture size number-pdf *n*(*a*) by
(24c)

Then, the mean trace length is given by: (24d)

Note that, for polydisperse fractures, .

Similarly for the horizontal pavement: (25a)

(25b)

and (25c)

If the outcrop window and pavement data are assumed to sample the same fracture population, the ratio between equations (25c) and (24d) yields: (26)

It is important to note that the result (equation 26) is valid whatever the fracture-size and orientation distributions.

Using equations (26) and (24b), the volumetric area is found as: (27)

where is the mean fracture area.

Then, the probability density function for the trace size is determined as: (28)

It should be noted that
is obtained as a numerical fit; thus, the mean fracture size <*a*> is determined from:
(29)

Finally, the measurement of *Σ*_{t,o} is used in equation (24a) to determine *ρ*:
(30)

where is found similarly to using the trace distribution on the pavement.

Results for all the data are given in Tables 2 & 3.

## Predictions of Percolation and Permeability

This section addresses the percolation and the permeability of the fracture networks by two different methods. The first method uses existing theoretical or semi-empirical solutions. The second one uses numerical calculations; numerical samples of fracture networks are generated according to the model parameters and the flow equations (equation 10) are solved.

### Snow’s formula

Since generic solutions exist only for a single Fisher family at the moment (Mourzenko *et al.* 2011*b*), the permeability of fracture networks that are described by at least two Fisher families can only be approximated by Snow’s formula (equation 12). This formula, which involves only the volumetric area of fractures and their orientation, regardless of their size and shape, becomes a good approximation when the density, *ρ*, is very large.

The geometrical tensor, ** Ψ** (equation 13b), can be directly obtained from the 3D orientation field data for each window, by expressing the normal vectors

**in a Cartesian referential (**

*n**x, y, z*) and averaging their dyadic product. For each window,

**is derived as: (31)**

*Ψ*where *α* is the angle between the fracture and the window; that is, the dyadic product of normals is corrected by the probability to observe the fracture.

For illustration, the fractures can be regarded as empty channels with a uniform aperture *b*_{f} = 1 mm. The transmissivity, σ, results from the cubic law:
(32)

Snow’s formula was applied to the total distribution (equation 16), and also applied separately to all of the windows that were analysed. Let us examine the results for windows 4 and 5 in Table 4. Three applications of Snow’s formula are listed in this table depending on the choice for the distribution. *f*_{tot} is given either by the measured distribution (equation 16), or by its approximations by two and four Fisher distributions. It is clear that these three estimations are not much different. The permeability for Window 4 is generally slightly smaller than for Window 5.

### Direct simulations

The measured parameters (Table 3) and hypotheses in the two previous sections are used first to generate reconstructed fracture networks that are statistically equivalent to the real ones.

The fractures are generated in a periodic cell of size *L _{x}×L_{y}×L_{z}*, with

*L*=60 m and

_{x}=L_{y}*L*60

_{z}=*/f*m, following our usual procedure (Huseby

*et al.*1997). The centres obey a Poisson distribution; that is without any spatial organization. The fractures are rectangles with the aspect ratio

*f*and their longer side is horizontal. The mean horizontal fracture size, <

*b*>, is kept constant for all the simulations, and the vertical one, <

*a*>, varies according to

*f*. Since the size distribution,

*n*(

*a*), depends on the parameter

*γ*, which is constant (equation 20), and on the mean fracture size <

*a*>, the aspect ratio is sufficient to determine

*n*(

*a*) for the numerical model. Then, the density,

*ρ*, which corresponds to each window, is used to reconstruct the fracture networks with a constant aperture, as discussed before.

However, the influence of the geometry of the fracture networks on the hydraulic properties is described by a universal parameter *ρ*′, the dimensionless density. Therefore, various values of *ρ*′ should be explored. This dimensionless density is the mean number of intersections per fracture. It is proportional to *ρ* and is calculated *a posteriori* on the generated networks, since equation (3a) applies only when the fractures are isotropically orientated within their plane.

The numerical simulations are performed with the following parameters, which cover the range of (*f, ρ*) measured in the outcrop windows (Table 3). The aspect ratio *f* is equal to 4, 5, 6, 7 and 12. For each *f*, the dimensionless density is *ρ*′=6.4, 7.4, 11.8, 15.6 and 23.6.

It should be noted that the network with (*f* =4, *ρ*′=11.8) is equivalent to Window 7, with (*f* =5, *ρ*′=15.6) is equivalent to Window 6, (*f* =6, *ρ*′=11.8) equivalent to Window 4, and (*f* =7, *ρ*′=15.6) to Window 5.

The fracture orientations are taken according to the probability distribution function *f*_{tot}(*θ, ψ*) (equation 16). However, the Fisher distributions given by two (equation 18) and four families (equation 17) are also used for a few cases.

In order to solve the flow problem, the fracture network is meshed by triangles. The resulting triangular mesh explicitly contains the fracture intersections. The mesh is constructed by an advancing front method, starting from the fracture contours and the intersection lines (Koudina *et al.* 1998). The prescribed spatial resolution is *δ*=0.25 m and corresponds to the maximum edge length of the triangles.

Then, the reconstructed networks are used to calculate the permeability tensor (Mourzenko *et al.* 2011*b*). The pressure is determined at the nodes of the mesh, and the flow equations are solved by use of a finite volume method with a first-order spatial discretization. The pressure gradient is regarded as uniform over each triangular surface element. The conservation equation (equation 10) is integrated over the control volumes around each mesh point. The resulting set of linear equations is solved by a conjugate gradient algorithm. Additional details for the numerical calculations of permeability are given by Adler *et al.* (2012).

A set of polydisperse fracture networks has been reconstructed as described previously. An example of the resulting 3D sample is shown in Figure 10a. The 60×60×8.6 m^{3} domain contains 14 603 fractures. This example corresponds to the fracture network equivalent to Window 5, with *f*=7 and *ρ*′=15.6. The *x*-axis points eastwards, the *y*-axis northwards and the *z*-axis upwards. The trace map obtained in a vertical plane orientated N140° as Window 5 is plotted in Figure 10b.

The flow equations are solved numerically on the local scale in the reconstructed polydisperse fracture networks. The permeability tensor, ** K**, is calculated. The mean of the diagonal components

*K*=tr(

**)/3 is shown in Figure 11 as a function of the aspect ratio,**

*K**f*, and the dimensionless density,

*ρ*′. It strongly increases as a function of

*ρ*′ and only slightly varies with

*f*. It should be noted that a few additional values of

*ρ*′=3.7, 4.3, 6.5, 8.5 and 13.8 are used for

*f*=12 to enable detailed reconstruction of the behaviour of

*K*.

The previous results were obtained for reconstructed fracture networks, with the fractures orientated according to the observed orientation probability distribution function *f*_{tot}(*θ, ψ*) (equation 16). For comparison, a few cases of networks reconstructed according to the model orientation distributions are considered below. The distribution modelled by two, *f*_{2F} (equation 18), and by four Fisher families, *f*_{4F} (equation 17), are tested for two cases that correspond to windows 4 and 5; that is, with (*f* =6, *ρ*′=11.8) and (*f* =7, *ρ*′=15.6), respectively. The results for windows 4 and 5 are given in Table 4.

It is seen that Snow’s formula provides, in all cases, a very good approximation to the overall results in the horizontal directions, although it clearly overestimates the vertical diagonal component of the permeability tensor. This deviation results from a greater tortuosity for vertical flow due to the horizontally elongated fracture shapes, which is unaccounted for in equation (12). The agreement improves when *ρ*′ increases or *f* decreases.

It should be emphasized that all of these estimations are proportional to the fracture transmissivity, *σ*, which was estimated by equation (32). This value has been used only to be able to provide the results in Darcies. However, Hamzehpour *et al.* (2009) showed that the mean field approximation provides an excellent estimation of the network permeability when the fracture transmissivity is variable; in other words, the average value of the fracture transmissivity can be used to calculate the network permeability with a high degree of precision. Finally, as stated previously, the numerical code can cope with a variable transmissivity at every node of the network.

## Conclusions

This paper presents methods that are applied to obtain the geometrical characteristics of fracture networks using readily available outcrop observations. The geometrical characteristics are used as model parameters in order to determine theoretical predictions of percolation and permeability by means of Snow’s formula, for the direct numerical simulation of fracture networks and for the estimation of permeability.

The analysis of the data available in different outcrop windows yields a unique description of the fracture-orientation distribution. This distribution can be modelled by two and four Fisher families, with good agreement with the measured 3D data. However, the distribution with four families yields the best approximation of the 3D observations. Then, the fracture-size distribution is obtained as a function of only one parameter, the aspect ratio.

Based on this analysis, fracture networks were reconstructed for various values of dimensionless density *ρ*′ and aspect ratio *f* according to the measured and modelled fracture-orientation distributions. They are all percolating since the investigated range 6.4 ≤*ρ*′≤ 23.6 is well above the percolation threshold.

Direct numerical simulations were performed to obtain the permeability tensor of the reconstructed fracture networks. The results were compared to the Snow predictions. Snow’s formula is an upper bound of *K*, with the minimal and maximal eigenvalues along the horizontal and vertical directions, respectively. Owing to the horizontal orientation of the longest fracture side, for *ρ*′/*f* ≤ 2, the minimal eigenvalue of *K* was found along the vertical direction and the maximal along the horizontal one. The results converge to the Snow prediction for larger *ρ*′.

A sensitivity analysis was conducted in order to examine the influence of the various parameters on the network properties. The permeability only slightly depends on the fracture aspect ratio, while the dimensionless density of fractures strongly influences the hydraulic behaviour of the fracture network.

## Acknowledgments

This work at UPCM and PPRIME was partly supported by a grant of ExxonMobil, which is gratefully acknowledged.

- © 2014 EAGE/The Geological Society of London