1 Introduction

Most rocks are fractured to one extent or another. Fractures that are more permeable than their host rock can act as preferential (or at least additional) pathways for fluid to flow through the rock, which is relevant in several areas of earth science and engineering, e.g. radioactive waste disposal in crystalline rock, exploitation of fractured hydrocarbon and geothermal reservoirs, or hydraulic fracturing (Bonnet et al. 2001; Neuman 2005; Salimzadeh and Khalili 2015; Tsang et al. 2015). In describing or predicting flow through fractured rock, the effective permeability of the rock, comprising rock matrix and a network of fractures, is a crucial parameter and may depend on several geometric properties of the fractures/networks, such as size, aperture, orientation, and fracture density.

It is possible to compute this effective permeability by numerically modelling flow through discrete fracture networks (DFN) and upscaling the results (e.g. Ahmed Elfeel and Geiger 2012; Lang et al. 2014). However, this numerical upscaling is computationally expensive. Considering the fact that the geometric information that is typically available on fracture networks is stochastic in nature, several realisations are required to obtain representative values for effective permeability. In addition, due to the degree of uncertainty inherent to fracture-network properties, one is often interested in quantifying the effect of this uncertainty on effective permeability by probing a multi-dimensional parameter space, e.g. when assessing the suitability of sites for radioactive waste disposal. Thus, there is a necessity for predictive analytical methods.

A prominent example of a predictive model in three dimensions is the heuristic equation proposed by Mourzenko et al. (2011). It was developed based on a vast collection of numerical simulations and theoretical arguments. Sævik et al. (2013) confirm the accuracy of the method, though it is, per definition, not applicable at low fracture densities when the background matrix has a non-negligible permeability.

The present work studies inclusion-based effective medium models which treat fractures as spheroidal inclusions embedded in a homogeneous and permeable rock matrix. Research on the use of such models to predict effective fractured-rock permeability, particularly at high fracture densities with intersecting fractures, has advanced rapidly recently (Fokker 2001; Pozdniakov and Tsang 2004; Barthélémy 2009; Sævik et al. 2013, 2014). Following the work of Sævik et al. (2013, 2014), the suitability and accuracy of some inclusion-based effective medium models are investigated here. In addition, the framework of the effective medium theory (EMT) is exploited to identify two characteristic regimes, depending on whether the overall flow through the rock mass is dominated by the fracture network or by the matrix. Due to distinctly differing flow patterns, the dependence of effective permeability on fracture density differs in the two regimes.

The effective permeability expressions obtained from the various effective medium theories (symmetric and asymmetric self-consistent, differential, Maxwell) are applied here to polydisperse fracture networks. Comparisons are also made with the Hashin–Shtrikman bounds, Snow’s equation, and Mourzenko’s equation.

2 Characterisation of Fractures and Fracture Networks

For the purpose of the effective medium models, a fracture is defined here geometrically as an oblate spheroidal inclusion within a rock matrix. The two equal semi-axis lengths are the fracture radius, R, and the length of the short semi-axis is one half of the maximum aperture, i.e. h / 2. Then, the aspect ratio, \(\alpha \), of the fracture is

$$\begin{aligned} \alpha = \dfrac{h}{2R} \quad (\alpha \ll 1) \end{aligned}$$
(1)

Several other conceptual models (e.g. Mourzenko et al. 2011), however, define fractures as flat objects (in this case, discs) with an equivalent aperture, \(a_{\mathrm {f}}\). These two geometrical representations can be related by equating the volumes of a single fracture, yielding:

$$\begin{aligned} a_{\mathrm {f}} = \frac{4}{3} R \alpha \end{aligned}$$
(2)

A fracture set is defined as a group of fractures which share common properties, e.g. radius and permeability. The density, \(\rho \), of fracture set i is the number of fracture centres, \(N_{\mathrm {f}}\), per volume, i.e.

$$\begin{aligned} \rho _i = \dfrac{N_{\mathrm {f},i}}{L^3} , \end{aligned}$$
(3)

where L is the length of a finite-sized cubic domain for which the permeability is to be determined. Note that, in this study, \(L>2R\). Let \(\phi _i\) be the total volume occupied by a fracture set per unit volume of rock mass. Assuming all fractures have the same aspect ratio,

$$\begin{aligned} \phi _i = \frac{4}{3} \pi \alpha R_i^3 \rho _i , \quad \phi = \sum _i^{{n}} \phi _i , \end{aligned}$$
(4)

where n is the total number of fracture sets. The following definition for fracture density is used:

$$\begin{aligned} \varepsilon _i = \rho _i R_i^3 , \quad \varepsilon = \sum _i^{{n}} \varepsilon _i \end{aligned}$$
(5)

Several other definitions of fracture density have been used in the literature (see, e.g. Mourzenko et al. 2011; Sævik et al. 2014; Li and Li 2015). The one used here is common, particularly in the context of effective medium models.

Each fracture is assigned a permeability, \(K_{\mathrm {f}}\), which may depend on several factors including aperture, surface roughness, filling material. The un-modified cubic law (Zimmerman and Bodvarsson 1996) is used here, for concreteness, and to provide a simple way for readers to relate the order of magnitude of the fractures permeability to its geometry:

$$\begin{aligned} K_{\mathrm {f}} = \dfrac{a_{\mathrm {f}}^2}{12} \end{aligned}$$
(6)

The cubic law is not a prerequisite for using these effective medium methods, and other expressions can be used for the permeability of the individual fractures, which attempt to account for factors such as roughness and contact areas (Zimmerman and Bodvarsson 1996). Moreover, a level of upscaling is inherent in the assignment of a “permeability” to the fracture, since the above equation is obtained by integrating the velocity across the entire aperture of the fracture, after which the “fracture” is treated as a locally homogeneous porous medium with a locally uniform permeability.

Finally, \(\kappa \) is the permeability ratio,

$$\begin{aligned} \kappa = \dfrac{K_{\mathrm {m}}}{ K_{\mathrm {f}} } \quad (\kappa \ll 1) , \end{aligned}$$
(7)

where \(K_{\mathrm {m}}\) is the permeability of the background rock matrix. In the context of this work, fracture permeability is always several orders of magnitude greater than the matrix permeability, i.e. \(\kappa \) is very small; the case of “filled” fractures having low permeability is not considered. The ratio, \(\alpha /\kappa \), between the two small, dimensionless parameters will be shown below to be an important characteristic of the system.

3 Upper Bounds

The Hashin–Shtrikman (H–S) bounds are rigorous upper and lower bounds for the effective conductivity of a two-component system. In the present context of fractured-rock permeability, the lower bound is almost equal to the matrix permeability, and is therefore trivial. The upper bound, however, yields interesting and non-trivial information. This upper bound is expressed here for monodisperse, isotropic networks as (Zimmerman 1989):

$$\begin{aligned} K^{+}_{\mathrm{HS}} = K_{\mathrm {f}} + \dfrac{3K_{\mathrm {f}} \left( K_{\mathrm {m}} - K_{\mathrm {f}}\right) (1-\phi )}{3 K_{\mathrm {f}} + \left( K_{\mathrm {m}} - K_{\mathrm {f}}\right) \phi } \end{aligned}$$
(8)

An early, simple model for the permeability of a fractured rock mass, which is only defined when \(K_{\mathrm {m}}=0\), was introduced by Snow (1969) under the assumption that the fractures in the domain are infinitely long and therefore are part of the percolating network. His estimate was extended by Oda (1995) to include fractures of finite length. This estimate will be referred to hereafter as Snow’s model. In the monodisperse, isotropic case, this estimate is

$$\begin{aligned} K_{\mathrm{Snow}} = \frac{2}{3} \rho \pi R^2 K_{\mathrm {f}} a_{\mathrm {f}} \end{aligned}$$
(9)

Using the definitions for equivalent aperture in Eq. (2) and fracture volume fraction in Eq. (4), Eq. (9) can be simplified to

$$\begin{aligned} K_{\mathrm{Snow}} = \frac{2}{3} \phi K_{\mathrm {f}} \end{aligned}$$
(10)

In order to compare this to the Hashin–Shtrikman upper bound, the H–S bound can be evaluated at zero matrix permeability, yielding

$$\begin{aligned} K^{+}_{\mathrm{HS}} = \dfrac{2\phi K_{\mathrm {f}}}{3-\phi } . \end{aligned}$$
(11)

Comparison of Eqs. (10) and (11) shows that

$$\begin{aligned} K^{+}_{\mathrm{HS}} = \dfrac{K_{\mathrm{Snow}}}{1-(\phi /3)} . \end{aligned}$$
(12)

Since the fracture porosity \(\phi \) is necessarily very small, it follows that Snow’s model essentially coincides with the H–S upper bound. This justifies the commonly made, but not previously proven, heuristic assumption that Snow’s model gives an “upper bound” to the actual permeability.

Fig. 1
figure 1

The parameter \(\beta \) for the dilute-limit solution [see Eq. (13)] plotted as a function of aspect ratio, \(\alpha \), for various values of the permeability ratio, \(\kappa \). Here, when \(\alpha \ll \kappa \), \(\beta \) does not depend on aspect ratio. Conversely, when \(\alpha \gg \kappa \), \(\beta \), and hence effective permeability, does not depend on fracture permeability. Though the latter statement is obviously not correct at high fracture densities with multiply intersecting fractures, the figure helps demarcate the two distinctly different regimes \(\alpha \ll \kappa \) and \(\alpha \gg \kappa \)

Full size image

4 Dilute-Limit Solution

The starting point for effective medium models is the calculation of the effect of a set of inclusions on the effective property at the dilute limit, i.e. an infinitesimally small concentration of inclusions (Zimmerman 1989). For randomly oriented spheroidal inclusions, this problem has been solved by Fricke (1924):

$$\begin{aligned} K&= K_{\mathrm {m}} \left( 1 + \beta \phi \right) \end{aligned}$$
(13)
$$\begin{aligned} \beta&= \dfrac{(r-1)}{3} \,\left[ \dfrac{4}{2 + (r-1)M} + \dfrac{1}{1 + (r-1)(1-M)} \right] \end{aligned}$$
(14)

were K is the effective permeability of the rock mass, \(r = \kappa ^{-1}\), and

$$\begin{aligned} M = \dfrac{2 {\varTheta } -\sin 2{\varTheta }}{2\tan {\varTheta } \sin ^2 {\varTheta }} , \quad \text {where } {\varTheta } = \arccos \alpha . \end{aligned}$$
(15)

For flat spheroids, i.e. when \(\alpha \ll 1\), \(M \approx \alpha \frac{\pi }{2}\) (Sævik et al. 2013, 2014) which is also a small number, leading to \(1-M \approx 1\). Combining these approximations with \(r~\gg ~1\), Eq. (14) becomes

$$\begin{aligned} \beta = \dfrac{4r}{3 \left( 2+\frac{\pi }{2}\alpha r \right) } = \dfrac{8}{3\kappa \left( 4+ {\pi }\frac{\alpha }{\kappa } \right) }. \end{aligned}$$
(16)

Figure 1 shows the variation of \(\beta \) with \(\alpha \) for various values of \(\kappa \). The graphs clearly show the two asymptotic regimes,

$$\begin{aligned} \beta = \left\{ \begin{array}{ll} \dfrac{8}{3 \pi \alpha } &{} \text {for } \alpha \gg \kappa \\ \\ \dfrac{2}{3 \kappa } &{} \text {for } \alpha \ll \kappa \end{array} \right. \end{aligned}$$
(17)

along with the nonlinear transition between the two. Note the apparent paradox that, when the aspect ratio is much greater than the permeability ratio, i.e. \(\alpha \gg \kappa \), the effective permeability does not depend on the fracture permeability. This peculiar fact alone highlights the need for effective medium models to extend the results to higher fracture densities, in which regime the effect of \(K_{\mathrm {f}}\) will eventually appear in the resulting expressions. Nonetheless, Fig. 1 already hints at the significance of \(\alpha /\kappa \) as a key parameter in determining effective permeability.

5 Effective Medium Models

The various effective medium models, as derived and presented by Sævik et al. (2013), are reviewed here. These include the asymmetric and symmetric self-consistent methods, the differential method, and Maxwell’s approximation. All of the methods are presented below for a finite number of fracture sets, n, and their various contributions are summed. However, for a continuous distribution of fracture radii, these contributions need to be integrated. In such cases, the integration can be performed numerically.

5.1 Maxwell Approximation

Following Maxwell’s conceptual model, the overall effect of several inclusions on the pressure distribution, and hence on permeability [Eq. (13), summed over all fracture sets], is equated to the effect of a large sphere of permeability K, representing the effective medium embedded in the same background rock matrix as the inclusions (see, e.g. Sevostianov 2014; Lutz and Zimmerman 2016). i.e.,

$$\begin{aligned} \sum _i^{{n}} \beta _i \phi _i = \beta _{\mathrm{sphere}} \end{aligned}$$
(18)

\(\beta _{\mathrm{sphere}}\) can be calculated from Eq. (14) with \(M_{\mathrm{sphere}}=2/3\) and \(r_{\mathrm{sphere}}=K/K_{\mathrm {m}}\). Rearranging the above equation yields

$$\begin{aligned} K = K_{\mathrm {m}} \dfrac{\left( 1+ \frac{2}{3}\sum _i^{{n}} \phi _i \beta _i \right) }{1 -\frac{1}{3} \sum _i^{{n}} \phi _i \beta _i} . \end{aligned}$$
(19)

Inserting the expression for \(\beta _i\) [Eq. (16)] for several fracture sets and \(\phi _i=\frac{4}{3} \pi \alpha \varepsilon _i\), one may rewrite the Maxwell equation as

$$\begin{aligned} K - K_{\mathrm {m}} = \dfrac{32}{9} \left( \dfrac{1}{3}K + \dfrac{2}{3} K_{\mathrm {m}} \right) \sum _i^{{n}} \dfrac{\varepsilon _i}{\dfrac{4 \kappa }{\pi \alpha } + 1} . \end{aligned}$$
(20)

5.2 Asymmetric Self-Consistent Method

Starting from the dilute-limit solution [Eq. (13) combined with Eq. (16)] for several fracture sets,

$$\begin{aligned} K -K_{\mathrm {m}}= \sum _i^{{n}} K_{\mathrm {m}} \dfrac{8}{3 \left( 4 \frac{K_{\mathrm {m}}}{K_{\mathrm {f},i}} + \pi \alpha \right) } \phi _i . \end{aligned}$$
(21)

This method attempts to account for the hydraulic influence of the various fractures on each other by replacing \(K_{\mathrm {m}}\) on the right-hand side of Eq. (21) by K. By also introducing fracture density, \(\varepsilon _i\), in the equation, one obtains

$$\begin{aligned} K -K_{\mathrm {m}}= \dfrac{32}{9} K \sum _i^{{n}} \dfrac{\varepsilon _i}{\dfrac{4 K}{\pi \alpha K_{\mathrm{f},{i}}} + 1} \end{aligned}$$
(22)

which is the form used in Sævik et al. (2013). Notice the implicit assumption that \(K_{\mathrm {f},i}\gg K\) in applying Eq. (16) here.

5.3 Symmetric Self-Consistent Method

This method is similar to the asymmetric self-consistent method, with the exception that the matrix is also represented as a set of inclusions, having spherical geometry (Sævik et al. 2013, 2014). Although this assumption seems arbitrary, and lacking a physical basis, it will be seen below that the predictions of this model are surprisingly accurate at low fracture densities. This model leads to the following expression (see Sævik et al. 2013):

$$\begin{aligned} K -K_{\mathrm {m}}= \dfrac{32}{9} \left( \frac{2}{3} K + \frac{1}{3} K_{\mathrm {m}}\right) \sum _i^{{n}} \dfrac{\varepsilon _i}{\dfrac{4 K}{\pi \alpha K_{\mathrm{f},{i}}} + 1} \end{aligned}$$
(23)

5.4 Differential Method

In this approach, fractures are assumed to be added in infinitesimal increments, with the result homogenised at each step, so that the \((j+1)\)th group of fractures are embedded in a hypothetical homogeneous medium whose permeability already reflects the effect of the first j fracture groups. Although the predictions of the differential method depend slightly on the integration path, for cases in which there are multiple inclusion sets, a natural integration path can be chosen such that the infinitesimal amount of each fracture set added, at each step, is proportional to the final proportion of the set’s volume fraction (Sævik et al. 2013), i.e.

$$\begin{aligned} (\phi _i/\phi )_j = (\phi _i/\phi )_{j+1} = \widehat{\phi _i}. \end{aligned}$$
(24)

The result is the following differential equation for the evolution of the effective permeability with fracture porosity (McLaughlin 1977; Zimmerman 1996; Sævik et al. 2013):

$$\begin{aligned} \dfrac{\text {d}K}{\text {d} \phi } = K \sum _i^{{n}} \beta _i {\widehat{\phi _i}} \end{aligned}$$
(25)

In the case of interest, where \(\alpha \) and \(\kappa \) are both very small, this equation takes the form (see Sævik et al. 2013)

$$\begin{aligned} \dfrac{\text {d}K}{\text {d} \varepsilon } = \dfrac{32}{9} K \sum _i^{{n}} \dfrac{{\widehat{\varepsilon _i}}}{\dfrac{4 K}{\pi \alpha K_{\mathrm{f},{i}}} + 1} , \end{aligned}$$
(26)

where \(\widehat{\varepsilon _i}\) is defined analogously to \(\widehat{\phi _i}\) in Eq. (24). For monodisperse, isotropic fracture networks, this differential equation can be integrated to yield (Sævik et al. 2013):

$$\begin{aligned} \frac{4}{\pi \alpha K_{\mathrm {f}}}( K - K_{\mathrm {m}} ) + \ln (K/K_{\mathrm {m}}) = \frac{32}{9} \varepsilon \end{aligned}$$
(27)

6 Heuristic Model of Mourzenko et al. (2011)

In addition to the effective medium methods, the equation for prediction of fracture-network permeability proposed by Mourzenko et al. (2011) (see also Bogdanov et al. 2003, 2007; Mourzenko et al. 2005) is included here for comparison. For polydisperse isotropic networks with disc-shaped fractures embedded in a permeable medium:

$$\begin{aligned} K&= K_{\mathrm {m}} + K_{\mathrm{Snow}}\times \nonumber \\&\quad \left( 1 - \dfrac{1}{1 + \frac{7}{3} \sigma '^{-0.7}} \left[ 1 - \dfrac{\widehat{\beta } (\rho ' - \rho '_{\mathrm{c}})^2}{\rho ' \left\{ 1 + \widehat{\beta } \left[ 1- (\rho ' - \rho '_{\mathrm{c}}) \right] \right\} } \right] \right) \end{aligned}$$
(28)

where \(K_{\mathrm{Snow}}=\frac{2}{3} \rho \pi \langle R^2 K_{\mathrm {f}} a_{\mathrm {f}} \rangle \), \(\rho '=\pi ^2 \rho \langle R^3 \rangle \), \(\rho '_{\mathrm{c}}=2.41\) represents the percolation threshold for disc-shaped fractures (Mourzenko et al. 2011), and \(\sigma '=\dfrac{K_{\mathrm{f},\mathrm{max}} \,a_{\mathrm{f},\mathrm{max}}}{R_{\mathrm{max}} K_{\mathrm {m{}_{\phantom {0_0}}}}}\) (with the subscript max denoting values for the fracture of maximum size). Finally, \(\widehat{\beta }\) is a fracture-shape-specific fitting parameter. To date, no \(\widehat{\beta }\) value for discs is available, and so the value for hexagons, \(\widehat{\beta }=0.18\), given by Mourzenko et al. (2011) is used here. Eq. (28) has been suggested for use only when \(\rho '\ge 4\).

For monodisperse networks, \(K_{\mathrm{Snow}}\) reduces to Eq. (9), \(\rho ' = \pi ^2 \varepsilon \), and

$$\begin{aligned} \sigma '=\dfrac{4}{3}\dfrac{\alpha }{\kappa } . \end{aligned}$$
(29)

7 Numerical Model and Simulations

The predictions of the various approximate models are now tested against explicit numerical simulations. These use the method presented by Lang et al. (2014) to calculate fractured-rock permeability as implemented in the modelling framework CSMP++. The approach works on unstructured finite-element meshes where fractures are represented explicitly as lower-dimensional objects, i.e. surfaces embedded in a volumetric domain which represents the host rock. The full permeability tensor of the fracture–matrix ensemble is obtained by means of volume averaging of the local pressure gradients and fluxes for three independent steady-state flow simulations, and the subsequent solution of an overdetermined system of equations. Calculating the permeability tensor over a restricted sub-volume (here, 90 % of the total volume) of the flow model allows the use of arbitrary fracture-network geometries that need not be restricted by periodicity or other constraints (for details, see Lang et al. 2014).

Each effective-permeability data point represents a statistically averaged value (here median) over 20–40 realisations. Of the three eigenvalues of the permeability tensor, the mean is calculated as a representative value. Figure 2 shows examples of realisations of fracture networks as used in the simulations to determine permeability.

Fig. 2
figure 2

Examples of varyingly dense polydisperse networks of disc-shaped fractures embedded in a cubic domain, as used in the simulations (from Case 5 in Table 1). The colouring shows fracture aperture which scales with radius in this example. As shown in (a), the element size used for the triangular mesh on the fractures depends on fracture size (see Table 1). a \(\varepsilon = 0.2\). b \(\varepsilon = 0.3\). c \(\varepsilon = 1.0\). d \(\varepsilon = 1.6\)

Full size image
Fig. 3
figure 3

Example of the dependence of the computed permeability \(K_\delta \) on spatial discretisation. Here, \(D=0.365\), \(R=10\,\mathrm{m}\), \(\varepsilon =1.12\), \(K=130.33 K_{\mathrm {m}}\)

Full size image

Each of these permeabilities can be corrected for discretisation errors in the same way as done by Mourzenko et al. (2011). A linear relationship between the computed permeability \(K_\delta \) and the dimensionless discretisation length \(\delta /R\) is assumed, where \(\delta \) is the maximum size of a triangular fracture element.

$$\begin{aligned} K = \dfrac{K_\delta }{1 + D \dfrac{\delta }{R}} \end{aligned}$$
(30)

For each set of simulations, D varies with fracture density and is determined at various densities by plotting the computed permeability as a function of \(\delta /R\) and extrapolating down to \(\delta /R = 0\) (see Fig. 3).

To explore the characteristics of the effective medium models, six different cases are defined here, each case posing a particular challenge to the prediction method. All of the fracture sets are isotropic and contained in a cubic domain with side lengths of 100 m. See Table 1 for a summary. Cases 1 and 2 are monodisperse networks. From the dilute-limit solution, it is clear that there is a substantial difference in the dependence of effective permeability on fracture density (or volume fraction) when the aspect ratio is much less or much greater than the permeability ratio. These two regimes, \(\alpha \ll \kappa \) and \(\alpha \gg \kappa \), are studied in Cases 1 and 2, respectively.

In Case 3, the fracture network comprises two isotropic fracture sets of strongly differing fracture radii (and, as such, aperture and permeability). Assuming that the relative number of fractures in each of the fracture sets follows a power-law distribution (many short fractures and few long fractures) with an exponent of \(n_{\mathrm{p}}=3.5\), the densities of the fracture sets are linked:

$$\begin{aligned} \rho _{2}=\rho _{1} \left( \frac{R_{1}}{R_{2}} \right) ^{3.5} \end{aligned}$$
(31)

Finally, Cases 4–6 consider continuous distributions of fracture radii, using the radii of the two fracture sets in Case 3 as the minimum and maximum values for a power-law distribution.

$$ \begin{aligned} P(R)&= \dfrac{1-n_{\mathrm{p}}}{R_{\mathrm{max}}^{1-n_{\mathrm{p}}} - R_{\mathrm{min}}^{1-n_{\mathrm{p}}}} R^{-n_{\mathrm{p}}}; \nonumber \\ \text {with }&n_{\mathrm{p}} = 3.5 \text { in Case 4}, \nonumber \\&n_{\mathrm{p}} = 1.5 \hbox { in Cases 5}\, \& \,6; \nonumber \\ \text { and }&4\,\mathrm{m} \le R \le 20\,\mathrm{m}. \end{aligned}$$

\(P(R)\text {d}R\) is the probability of observing fracture radii within \([R,R+\text {d}R]\) (see, e.g. Lang et al. 2014; Mourzenko et al. 2005).

8 Comparison to Numerical Data

A comparison between the various effective medium models, the upper bound(s), Mourzenko’s model, and numerical data is achieved here by plotting \(K/K_{\mathrm {m}}\) as a function of \(\varepsilon \). For the numerical data, the range between the 5th and 95th percentiles of the 20–40 realisations are shown as a measure of spread in addition to the median value (see Fig. 4).

8.1 Case 1

As shown in Fig. 4a, all of the effective medium methods, as well as Mourzenko’s equation, give accurate predictions for effective permeability in this regime, except the Maxwell and symmetric self-consistent methods. Their predictions start to deviate from the data when \(\varepsilon >0.5\), each in a different direction. The H–S/Snow upper bounds are very close to the actual permeability, which means that almost all the fractures within the medium contribute to flow.

Table 1 Simulation parameters for Cases 1–6
Full size table
Fig. 4
figure 4

Effective fractured-rock permeability: asymmetric self-consistent (SSC), asymmetric self-consistent (ASC), differential (Diff.), Maxwell (Maxw.), Mourzenko et al. (2011) (Mourz.), H–S and/or Snow’s upper bounds (HS+, Snow), and numerical data (Data). The upper bounds in c)–f) are from Snow’s model. The numerical data—represented by the median over 20–40 realisations, and the range between the 5th and 95th percentiles—comprises a total of 1431 realisations (mean permeabilities averaged over three eigenvalues), yielding 53 data points (median permeabilities). The maximum number of fractures in the results shown here range from 1375 (Case 5/6) to 10,035 (Case 3). a Case 1 (monodisperse, \(\alpha /\kappa =0.078\)). b Case 2 (monodisperse, \(\alpha /\kappa = 78\)). c Case 3 (two fracture sets, \(\left<\alpha /\kappa \right>_\phi =105\)). d Case 4 (power law with \(n_{\mathrm{p}}=3.5\), \(\left<\alpha /\kappa \right>_\phi =111\)). e Case 5 (power law with \(n_{\mathrm{p}}=1.5\), \(\left<\alpha /\kappa \right>_\phi =177\)). f Case 6 (power law with \(n_{\mathrm{p}}=1.5\), \(\left<\alpha /\kappa \right>_\phi =0.177\)) \(\left<\alpha /\kappa \right>_\phi \) is a volume-weighted average of \(\alpha /\kappa \) over the fracture sets (explained in Sect. 9.3.2)

Full size image

Neither the numerical data nor the approximate models predict a distinct percolation threshold here. Apart from Maxwell’s approximation and the symmetric self-consistent method, the results even appear to be linear. For this monodisperse case, exploiting \(\alpha /\kappa \ll 1\) in the asymmetric self-consistent [Eq. (22)] and differential [Eq. (27)] methods yields a linear relationship between K and \(\varepsilon \), irrespective of the value of \(\varepsilon \):

$$\begin{aligned} K = K_{\mathrm {m}} + \dfrac{8 \pi }{9} \alpha K_{\mathrm {f}} \varepsilon \end{aligned}$$
(32)

The second term on the right side of the above equation represents the effective permeability of the fracture network, if the matrix permeability were zero. Hence, the above result shows that the effective permeability of the fractured/porous rock mass is simply equal to the sum of the permeability of the matrix, plus the permeability of the “fracture network”. This type of simple additivity is very rare in effective medium problems and arises here only in the extreme range of the parameter space in which \(\alpha \ll 1\) and \(\kappa \ll 1\).

8.2 Case 2

Here, the regime \(\alpha \gg \kappa \) is studied. As in the opposite regime (Case 1), Fig. 4b shows that the H–S upper bound and the Snow model are practically equal. The Maxwell approach tends to infinity when the denominator in Eq. (19) tends to zero, i.e. when \(\varepsilon \rightarrow \frac{27}{32} \left( \frac{4 \kappa }{\pi \alpha }+1\right) \). Hence, in this particular case, Maxwell’s equation gives meaningless predictions when \(\varepsilon \ge 0.86\). In theory, the percolation threshold for this set of randomly oriented disc-shaped fractures is at \(\varepsilon =0.244\) (Mourzenko et al. 2005). As mentioned in Sævik et al. (2013), the asymmetric self-consistent method is able to predict the percolation threshold. However, the symmetric self-consistent approach is better at capturing the general behaviour of the curve at low fracture densities. At higher densities (\(\varepsilon >0.5\)), the same type of deviation from the numerical data as in Case 1 can be observed. In this range of higher densities, the gradient of the asymmetric self-consistent curve is more accurate, though it overestimates the effective permeability. The predictions of the differential method are much lower than the numerical data. Sævik et al. (2013) state that the differential method only gives useful results when \(\alpha \le \kappa \).

In this example, Mourzenko’s equation provides an excellent match to the numerical data, accounting well for both medium- and high-fracture-density ranges. Mourzenko et al. (2011) point out that, according to percolation theory, K is expected to be a quadratic function of \(\varepsilon \) at densities above but close to the percolation threshold. This statement is valid when \(K_{\mathrm {m}} \approx 0\), but also appears to hold here. At higher densities, K then becomes a linear function of \(\varepsilon \).

It can easily be shown that both self-consistent methods become linear in \(\varepsilon \) when \(K\gg K_{\mathrm {m}}\):

$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&\propto \dfrac{\pi }{4} \dfrac{\alpha }{\kappa } \left( \dfrac{32}{9}\varepsilon - 1 \right) \quad \text {asymmetric} \end{aligned}$$
(33)
$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&\propto \dfrac{\pi }{4} \dfrac{\alpha }{\kappa } \left( \dfrac{64}{27}\varepsilon - 1 \right) \quad \text {symmetric} \end{aligned}$$
(34)

Under the assumption that, at very low densities (i.e. below the percolation threshold), \(K \ll \alpha K_{\mathrm {f}}\) (which, in this case, implies that effective permeability is similar in magnitude to matrix permeability), one can derive the following equations from the effective medium models:

$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&=\dfrac{1}{1-\frac{32}{9} \varepsilon } \quad \text {asymmetric self-consistent} \end{aligned}$$
(35)
$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&=\dfrac{1+\frac{32}{27} \varepsilon }{1-\frac{64}{27} \varepsilon } \quad \text {symmetric self-consistent} \end{aligned}$$
(36)
$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&=\exp \left( \frac{32}{9} \varepsilon \right) \quad \text {differential} \end{aligned}$$
(37)
$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&=\dfrac{1+\frac{64}{27} \varepsilon }{1-\frac{32}{27} \varepsilon } \quad \text {Maxwell} \end{aligned}$$
(38)
$$\begin{aligned} \dfrac{K}{K_{\mathrm {m}}}&=1 + \frac{32}{9} \varepsilon \quad \text {dilute limit} \end{aligned}$$
(39)

In these expressions, the effective permeability is independent of fracture permeability because, well below the percolation threshold, the fluid always has to pass through the matrix, making matrix permeability the limiting factor.

8.3 Case 3

There are two fracture sets with distinctly different fracture radii (4 and 20 m) and, hence, differing fracture permeabilities. This produces a significant skew in the effective-permeability distribution, as can be seen in the asymmetric spread of the values of the numerical data in Fig. 4c. While the fracture set with many small fractures leads to a relatively low median value, the large-radius fracture set is responsible for the upper outliers. In general terms, the accuracy of predictions of the effective medium models does not seem to change because of this. They tend to over- or underestimate effective permeability in much the same way they do in Case 2. However, the skewed nature of the fracture network appears to pose a problem for Mourzenko’s heuristic model, which overestimates effective permeability in a manner very similar to the asymmetric self-consistent method.

8.4 Cases 4–6

The comparison between the effective medium models and numerical data using a power-law distribution (Fig. 4d–f) reinforces the findings of Cases 1–3. The symmetric self-consistent model performs best at low fracture densities and eventually diverges from the data. The asymmetric self-consistent model generally overestimates effective permeability, but converges towards the correct gradient of K to \(\varepsilon \) at large fracture densities (this is more obvious in Fig. 4e than in Fig. 4d), where full convergence has not yet occurred. The accuracy of both self-consistent methods is unaffected by the properties of the fracture-size distribution.

Mourzenko’s heuristic model performs very well in these cases. Its accuracy appears to increase with decreasing values of the power-law exponent, \(n_{\mathrm{p}}\). This is because part of the fitting of Eq. (28) relies on the dimensionless fracture transmissivity \(\sigma '\) of the largest fracture. As \(n_{\mathrm{p}}\) decreases, the influence and relevance of \(\sigma '\), and hence the accuracy of the fit, increases.

9 Discussion

Figure 5 summarises the various characteristics of effective rock-mass permeability addressed in this study. The ratio \(\alpha /\kappa \) is a measure of the relative contributions of the fracture network and the matrix to overall flow, going from no preferential flow through the fractures when \(\alpha /\kappa \rightarrow 0\) to no flow through the matrix when \(\alpha /\kappa \rightarrow \infty \).

Fig. 5
figure 5

Overview of flow characteristics and effective permeability of a fractured-rock mass. Figures depict slices through the middle of a 3D cubic domain from Cases 5 and 6 (right and left column, respectively). Contour lines and map represent pressure with flow from left to right. Fracture trace lines are shown as thick black lines. The given equations are derived from the asymmetric self-consistent method for monodisperse, isotropic networks

Full size image

In the following, one of the two self-consistent methods [Eq. (22) or (23)] is used to explore particular aspects of the permeability of a fractured rock mass.

9.1 Percolation Threshold

While there is a clear percolation threshold in Fig. 4b where \(\alpha \gg \kappa \), it is clearly lacking in Fig. 4a where \(\alpha \ll \kappa \). A closer look at this behaviour for various values of \(\alpha /\kappa \) and \(\varepsilon \le 0.5\) using the symmetric self-consistent method reveals how the linear behaviour of effective permeability with \(\varepsilon \) for \(\alpha /\kappa \ll 1\) (described in Sect. 8.1) transitions into the distinctly nonlinear regime as \(\alpha /\kappa \) increases (Fig. 6). Hence, the smaller the contribution of the rock matrix to flow, the more distinct the percolation threshold is. This occurs because, when the permeability of the matrix is non-negligible, fluid can easily bridge any gaps between unconnected fractures by passing through the matrix (Bogdanov et al. 2003).

Fig. 6
figure 6

Effective permeability at low fracture densities for various values of \(\alpha /\kappa \) using the symmetric self-consistent method. The parameters used here are based on Case 2 (see Table 1). Variations are achieved by varying the fracture radius, and correspondingly, fracture aperture and permeability. The two vertical lines represent the percolation threshold when the matrix is impermeable as determined numerically by Mourzenko et al. (2005) (dashed) and predicted by the symmetric self-consistent method (dash-dotted)

Full size image
Fig. 7
figure 7

Contour map showing the relative difference between the predictions of the asymmetric self-consistent method, \(K_\text {ASC}\), and the H–S upper bounds, i.e. \((K_{\mathrm{HS}}^{+} - K_\text {ASC})/K_{\mathrm{HS}}^{+}\). The parameters used here are the same as in Fig. 6

Full size image

In the numerical data generated by Bogdanov et al. (2007) of effective permeability for fracture networks and varying values of \(\alpha /\kappa \), the percolation threshold as a distinct, discontinuous feature tends to disappear when \(\alpha /\kappa \le 7.5\). This is in agreement with the findings discussed here, although Bogdanov et al. (2007) do not consider the case \(\alpha /\kappa \ll 1\).

9.2 Upper Bounds

Figure 7 is a contour plot of the relative difference between the H–S upper bounds and the effective permeability as calculated using the asymmetric self-consistent method for various values of \(\varepsilon \) and \(\alpha /\kappa \). It can be seen as a measure of the fraction of fractures that contribute to flow. When \(\alpha \gg \kappa \), this difference is maximal below the percolation threshold and reduces to a constant value with increasing \(\varepsilon \). Essentially, the contribution of non-percolating clusters to overall flow is very low. At very high densities, almost all clusters are connected, and only peripheral parts of the fractures do not contribute to flow (Mourzenko et al. 2011; Leung and Zimmerman 2012).

For \(\alpha \ll \kappa \), effective permeability is consistently close to the upper bounds. Flow through the matrix is important, always connecting fracture clusters and maximising the number of fractures that contribute to flow. This suggests that, in this regime, theoretical upper bounds are reasonably accurate estimates for effective permeability.

9.3 Polydisperse Networks

9.3.1 Smallest Fracture Radius

Polydisperse fracture networks are typically dominated by large fractures with high permeabilities. However, they also usually contain many small fractures with low permeabilities. A threshold fracture radius, \(R_{\mathrm{th}}\), may exist, for which, the contribution of smaller fractures to flow is negligible (see, e.g. Berkowitz et al. 2000; Dreuzy et al. 2001). When characterising fracture networks, it is important to know which fracture sets may be neglected within a given range of accuracy. The additive nature of effective medium methods provides a convenient framework for addressing this issue.

Consider, as an example, Case 4 in Sect. 8.4. Assuming now that there is no clear minimum fracture radius, i.e. \(R_{\mathrm{min}} = 0\,\mathrm{m}\) (as opposed to 4 m in Case 4), one may be interested in determining the threshold radius, \(R_{\mathrm{th}}\). This can be assessed, e.g. using the asymmetric self-consistent method, by neglecting the contributions of all fractures in Eq. (22) for which \(R_i < R_{\mathrm{th}}\). Figure 8 shows the effect this would have on the predicted effective permeability for various values of \(R_{\mathrm{th}}\). Note that \(\varepsilon \) in Fig. 8 (i.e. the x-axis) is calculated with Eq. (5) for all fracture radii.

Fig. 8
figure 8

Effect of neglecting the contribution of fractures with radii smaller than \(R_{\mathrm{th}}\) on effective permeability as calculated with the asymmetric self-consistent method. This is equivalent to Case 4 with \(R_{\mathrm{min}}=0\,\mathrm{m}\)

Full size image

9.3.2 Characteristic \(\langle \alpha /\kappa \rangle \)

For fracture networks with several fracture sets of differing \(\alpha /\kappa \), as is the case in polydisperse networks (including Cases 3–6), an obvious question that arises is how to define an average \(\langle \alpha /\kappa \rangle \) that is representative of the fracture network and rock matrix. A simple arithmetic average does not make sense, since it lacks information on the sizes of the fractures. Instead, a volume-weighted average seems more appropriate:

$$\begin{aligned} \langle \dfrac{\alpha }{\kappa } \rangle _\phi = \dfrac{1}{\phi }{\sum \limits _i^{\mathrm{n}} \dfrac{\alpha _i}{\kappa _i} \phi _i} \end{aligned}$$
(40)

Using the example in Sect. 9.3.1 (i.e. Case 4 with \(R_{\mathrm{min}}=0\,\mathrm{m}\)), effective permeability is plotted in Fig. 9 over fracture density. By varying \(R_{\mathrm{max}}\) between 0 and 20 m, various curves can be plotted covering a range of values of \(\langle \alpha /\kappa \rangle _\phi \). Here, the asymmetric self-consistent method is used to calculate K.

Fig. 9
figure 9

Effective permeability as a function of fracture density for various values of \(\langle \alpha /\kappa \rangle _\phi \) using the asymmetric self-consistent method. The parameters used here are based on Case 4 (see Table 1), with \(R_{\mathrm{min}}=0\,\mathrm{m}\). Here, \(R_{\mathrm{max}}\) is varied between 0 and 20 m

Full size image

Equation (40) can be used to calculate \(\langle \alpha /\kappa \rangle _\phi \) in Cases 3–6. In Case 3, this yields \(\langle \alpha /\kappa \rangle _\phi = 105.\) Cases 4–6 require an integration over the fracture radii, eventually leading to

$$\begin{aligned} \langle \dfrac{\alpha }{\kappa } \rangle _\phi = \dfrac{4 \alpha ^3 (4-n_{\mathrm{p}}) \left[ R_{\mathrm{max}}^{(6-n_{\mathrm{p}})} - R_{\mathrm{min}}^{(6-n_{\mathrm{p}})}\right] }{27 K_{\mathrm {m}} (6-n_{\mathrm{p}})\left[ R_{\mathrm{max}}^{(4-n_{\mathrm{p}})} - R_{\mathrm{min}}^{(4-n_{\mathrm{p}})}\right] } , \end{aligned}$$
(41)

i.e. \(\langle {\alpha /\kappa } \rangle _\phi = 111\) in Case 4, and \(\langle {\alpha /\kappa } \rangle _\phi = 177\) in Case 5. Finally, \(\langle {\alpha /\kappa } \rangle _\phi = 0.177\) in Case 6. These values adequately depict the nature of the curves shown in Fig. 4c–f, i.e. linear when \(\langle {\alpha /\kappa } \rangle _\phi = 0.177\) and very nonlinear for \(\langle {\alpha /\kappa } \rangle _\phi = 111\) and \(\langle {\alpha /\kappa } \rangle _\phi = 177\).

Figure 10 compares effective permeability predictions by the asymmetric self-consistent method for various fracture-size distributions. Here, a monodisperse case (Case 2) with \({\alpha /\kappa }=78\) is plotted on the same graph as two polydisperse cases (Cases 4 and 5). In Cases 4 and 5, the aspect ratio, \(\alpha \), has been altered such that \(\langle {\alpha /\kappa } \rangle _\phi =78\). The comparison touches on some of the points that have previously been made for the asymmetric self-consistent method when \(\alpha /\kappa >1\):

  1. 1.

    Below the percolation threshold, K does not depend on \(\alpha /\kappa \) [see Eq. (35)].

  2. 2.

    For a given value of \(\alpha /\kappa \): When \(K \gg K_{\mathrm {m}}\), K is linearly proportional to \(\varepsilon \) [see Eq. (33)].

However, one can also see that while \(\alpha /\kappa \) and \(\langle {\alpha /\kappa } \rangle _\phi \) can be used to characterise the system, they cannot be expected to uniquely represent the behaviour of the effective permeability, particularly when \(\alpha /\kappa >1\).

Fig. 10
figure 10

Effective permeability as a function of fracture density for various fracture-size distributions using the asymmetric self-consistent method. The three curves represent Case 2 (\({\alpha /\kappa }=78\)), and Cases 4 and 5 with altered values of \(\alpha \) to ensure that \(\langle {\alpha /\kappa } \rangle _\phi =78\) in both cases

Full size image

10 Previous Work Dealing with the Parameter \(\alpha /\kappa \)

Bogdanov et al. (2003) introduced a parameter to characterise their numerical data on monodisperse, isotropic networks,

$$\begin{aligned} \sigma '=\dfrac{K_{\mathrm {f}} \,a_{\mathrm {f}}}{R K_{\mathrm {m}}} , \quad \text {i.e.} \quad \sigma '=\dfrac{4}{3}\dfrac{\alpha }{\kappa } \end{aligned}$$
(42)

in the current notation. Their values for \(\sigma '\) ranged from \(10^{-6}\) to \(10^4\). However, in the cases where \(\sigma '~<~1\), fracture permeabilities were either similar in magnitude to or lower than matrix permeability, i.e. \(\kappa \approx 1\) or \(\kappa >1\). In the present work, however, only fractures that are very permeable compared to the matrix are considered, i.e. \(\kappa \ll 1\), making a comparison in the regime \(\sigma '<1\) infeasible.

Nonetheless, when Bogdanov et al. (2003) gradually increase \(\sigma '\) from 10 to \(10^4\), they also find an increasingly more distinct percolation threshold (as discussed in Sect. 9.1). Their investigation is limited to fracture densities of \(\varepsilon <1\). Bogdanov et al. (2003) seem to assume that the effective permeability is only linear with fracture density at very low fracture densities below the percolation threshold or when \(\kappa >1\). As pointed out in Sect. 8.1, if \(\alpha \ll \kappa \), this linearity also holds at high fracture densities even when \(\kappa \ll 1\).

In Bogdanov et al. (2007), the range \(1\le \sigma '\le 10^4\) is studied, again for \(\varepsilon <1\). Note that in the calibration of the model by Mourzenko et al. (2011), higher fracture densities are considered. Bogdanov et al. (2007) extend the definition of \(\sigma '\) to account for polydisperse networks and size-dependent fracture permeability (see Sect. 6):

$$\begin{aligned} \sigma '=\dfrac{K_{\mathrm{f},\mathrm{max}} \,a_{\mathrm{f},\mathrm{max}}}{R_{\mathrm{max}} K_{\mathrm {m}}} \end{aligned}$$
(43)

This depends solely on the properties of the largest fracture and is inherently different from the averaging conducted in Sect. 9.3.2.

Sævik et al. (2013) use the parameter

$$\begin{aligned} \lambda ^{-1} = \dfrac{\alpha }{\kappa } K_{\mathrm {m}} \end{aligned}$$
(44)

in their study, and varied \(\alpha /\kappa \) between 1 and 100. Here, again, \(\varepsilon <1\) (when \(K_{\mathrm {m}} \ne 0\)). Note that they considered both isotropic and anisotropic (albeit monodisperse) fracture distributions. Their results show a similar transition from an almost linear, to an increasingly nonlinear, character as depicted in Fig. 6.

Sævik et al. (2014) consider the parameter \(\lambda ^{-1} = \alpha /\kappa \), varying it between 10 and \(10^4\), with \(\varepsilon < 5\).

11 Conclusions

In predicting the effective permeability of a three-dimensional fractured rock mass, two distinct regimes can be distinguished, depending on the relative size of two small, dimensionless numbers: the aspect ratio of the spheroidal fractures, \(\alpha \), and the permeability ratio, \(\kappa \). When \(\alpha /\kappa \ll 1\), effective permeability is linearly dependent on fracture density without a distinct percolation threshold. With increasing \(\alpha /\kappa \), this relationship becomes increasingly nonlinear at low and intermediate fracture densities, but remains linear at high fracture densities. The ratio \(\alpha /\kappa \) can therefore be interpreted as a measure of the relative contributions of the fracture network and matrix to the overall flow, similar to the flux ratio discussed by Paluszny and Matthai (2010), and Nick et al. (2011). For polydisperse fracture networks, a characteristic value of \(\langle \alpha /\kappa \rangle \) can be determined by weighting the individual values of \(\alpha /\kappa \) of the various fracture sets with fracture volume while averaging.

Comparison to explicit numerical simulations of mono- and polydisperse isotropic fracture networks (with fracture-size dependent fracture permeabilities) shows that the self-consistent effective medium methods are generally capable of predicting effective permeability over a wide range of conditions. While the symmetric self-consistent method is particularly accurate at low fracture densities, the asymmetric self-consistent method predicts the correct asymptotic behaviour (cf. Sævik et al. 2013). The differential method is only useful when \(\alpha /\kappa <1\). In that regime, even the Hashin–Shtrikman and Snow upper bounds appear to give good approximations. Maxwell’s approximation is only reliable at very low fracture densities.

These effective medium models have been shown to be powerful tools. They are valid at both low and high fracture densities. When considering polydisperse fracture networks, the type and characteristics of the fracture-size (as well as aperture and fracture-permeability) distribution do not affect the applicability or accuracy of these models. These results are perhaps surprising, since intuitively, fracture intersections would seem to be a crucial factor in controlling the effective permeability, and yet these models do not explicitly contain fracture intersections as a parameter. Moreover, they are all based on the problem of a single “inclusion”, for which the concept of intersection is not even meaningful.

It should be noted that the heuristic model proposed by Mourzenko et al. (2011) is good at predicting effective permeability. It is very accurate for monodisperse fracture networks, but suffers a slight loss of accuracy for polydisperse networks, that appears to depend on the attributes of fracture-size distribution. The present work has in effect shown a link between predictions of methods based on inclusions-in-a-matrix and methods based on fracture networks.

Finally, for the case of zero matrix permeability, the well-known approximation of Snow, which is based on network considerations rather than a continuum approach, is shown to essentially coincide with the upper Hashin–Shtrikman bound, thereby proving that the commonly made assumption that Snow’s equation is an “upper bound” is indeed correct.