On the Density Variability of Poissonian Discrete Fracture Networks, with application to power-law fracture size distributions

. This paper presents analytical solutions to estimate at any scale the fracture density variability associated to stochastic Discrete Fracture Networks. These analytical solutions are based upon the assumption that each fracture in the network is an independent event. Analytical solutions are developed for any kind of fracture density indicators. Those analytical solutions are veriﬁed by numerical computing of the fracture density variability in three-dimensional stochastic Discrete Fracture Network (DFN) models following various orientation and size distributions, including the heavy-tailed power-law fracture size distribution. We show that this variability is dependent on the fracture size distribution and the measurement scale, but not on the orientation distribution. We also show that for networks following power-law size distribution, the scaling of the three-dimensional fracture density variability clearly depends on the power-law exponent.


Introduction
Characterizing fracture networks in geosciences is a key challenge for many industrial projects such as deep waste disposal, hydrogeology or petroleum resources, because it may change the mechanical (Davy et al., 2018;Grechka and Kachanov, 2006) and hydrological (Bogdanov et al., 2007;De Dreuzy et al., 2001a, b) behaviour of the rock mass. Fractures being ubiquitous and at all scales, the description of these physical properties is often far beyond the reach of a continuum approach (Jing, 2003). Discrete Fracture Networks are computational models explicitly representing the geometry of fractures in a network, and can be used as a basis for physical simulations (mechanical strength, flow, transport. . . ) (see Jing, 2003;Lei et al., 2017, for a review). Considering the scarce nature of geological data, statistical methods have been widely used to generate DFN models, where all fracture geometrical attributes are treated as independent variables from probability distribution derived from the field. Indeed, fracture networks are often described from size distribution, sets of orientations, location, and densities (Dershowitz and Einstein, 1988). Unfortunately, the difficulty of access and resolution to volumetric data makes it difficult to directly measure three-dimensional fracture densities. Stereological analysis proposes theoretical relationship to calculate the 3-D density from 1-D or 2-D measurements under some assumptions (Berkowitz and Adler, 1998;Darcel et al., 2003a;Warburton, 1980). For example, the total fracture surface per unit volume p 32 can be calculated from one-dimensional fracture intersections along scanlines or boreholes using some conversion factor (Mauldon, 1994;Terzaghi, 1965;Wang, 2005). Most of these studies focus on characterizing the mean density of a fractured system with poor interest to the underlying variability. However, fracture density variability is an indicator for fracture clustering, which may have dramatical impact on connectivity (Darcel et al., 2003b;La Pointe, 1988;Manzocchi, 2002). Darcel et al. (2013) proposed an analytical solution to quantify at any scale, the standard deviation associated to fracture frequency p 10 , considering fracture-borehole crossing as a one-dimensional Poisson point process (random positions with fixed density). They show that the standard deviation associated to the number of intersections per unit length p 10 is inversely proportional to the square root of the measurement scale. This solution was also demonstrated by Lu et al. (2017), who validated it numerically computing p 10 values on boreholes crossing three-dimensional Poissonian Published by Copernicus Publications on behalf of the European Geosciences Union.
DFN models of fixed fracture density. In this paper, we aim to develop analytical solutions for three-dimensional Poissonian DFN models to quantify the range of variability of any kind of fracture density, of any dimension. Particularly, we propose solutions for three-dimensional density variability that cannot be obtained directly from the field. We show that this variability depends on the measurement scale and the fracture size distribution, but not on the orientation distribution. These solutions are validated by numerical simulations, computing the associated fracture densities mean and variance at various scales for three-dimensional Poissonian Discrete Fracture Networks.

Theoretical development
In the following, we will calculate the variability of the fracture densities p Dd such as defined by Dershowitz and Herda (1992) and Mauldon and Dershowitz (2000). D and d refers to the dimensions of the embedding system and of the fracture trace in , respectively, so that d ≤ D. For instance, p 10 is the number of intersection points (d = 0) along a line (D = 1) per unit line size. We denote by µ f the contribution of a fracture f to p Dd (µ f is a d-dimension measure) and s the scale at which p Dd and µ f are calculated so that: First, we calculate the spatial variability of the measure µ f (s) calculated in a D-sphere of size s embedded in a larger space of volume V . In D-dimension, a fracture is represented by an object of dimension D −1 (a point in a scanline, a line in a 2-D outcrop, or a surface in the 3-D space) with a size S f . The proportion of the total volume occupied by the fracture -i.e., the volume where the measure µ f (s) is not nil -is: We then assume that, when it is not nil, the measure µ f (s) derives from a distribution with a mean µ f,s and a standard deviation σ µ f,s . For the total system, we can calculate the first and second order moment of µ f (s): The variance σ 2 (µ f ) is written as: We now sum the contribution of all fractures to calculate µ (s) = f µ f (s). With the Poisson's hypothesis, all fractures are independent events and the variance σ 2 (µ) is the sum of the variance of each fracture: And finally, the variability of the density measure p Dd is σ 2 (p Dd ) = σ 2 (µ)/s 2D . In the following, we compute the dimensionless ratio between the variance σ 2 (p Dd ) and the square of the mean p Dd : In the fractal theory, λ (·) is called lacunarity (Plotnick et al., 1996). This is a scale dependent measure, whose analysis gives an idea of the scaling of fracture network textural heterogeneity and potentially on different regimes (Plotnick et al., 1996;Roy et al., 2014). For geological environments, fracture networks are characterized by a wide distribution of fracture sizes. We denote by n (l) the density distribution of fracture sizes, and n (l) dl is the number of fractures of size in the range [l, l + dl] per system volume (Davy, 1993). The total number of fractures in a system of volume V is N = V n (l) dl . We then assume that the measure µ f (s) and the occupied volume ratio P f (s) only depends on the fracture size l so that µ f (s) = µ(l, s) and P f (s) = P (l, s). Equation (6) can now be written as a function of the fracture size distribution: In the following, we especially focus on the power-law size distribution, which have been found to adequately fit natural systems (Bonnet et al., 2001;Bour, 2002;Bour and Davy, 1999;Odling, 1997): with a the power-law exponent, α the density term.

1-D measurements
Fracture abundance is often quantified on boreholes using the one-dimensional fracture frequency p 10 measurement, defined as the number N of crossing fractures per unit borehole length L. If a fracture intersecting this borehole is also part of a subsample of size s, the associated measure µ(l, s) = 1. The ratio of subsamples intersected by this fracture is P (l, s) = s/L. Considering that s L, then P (l, s) 1, and the lacunarity of such measurement is the one of a onedimensional Poisson point process (Darcel et al., 2013;Lu et al., 2017): + σ 2 ( (l, s)) and P (l, s).

3-D measurements
In three-dimensional systems, the fractures d-measures µ f are either their occurrence (d = 0), surface (d = 2) or surrounded sphere (d = 3). µ f (s) such as defined in the previous section is the part of the measure µ f that is included into a domain S of size s, that is the product of µ f by the occurrence probability of the fracture to be included into the domain S, f,s (also noted (l, s) following the notation of the previous section). We define f,s as the ratio between the fracture surface included into the domain S to the total fracture surface. If the fracture size is larger than s, both the average value of f,s and its standard deviation are proportional to the surface s 2 . If P (l, s) 1, that is for large systems, Eqs. (4) and (7) requires to know the expression of (l, s) 2 + σ 2 ( (l, s)). The scaling of (l, s) 2 + σ 2 ( (l, s)) and P (l, s) are given in Table 1, with β s and β f shape factors depending on domain S and fractures geometries respectively. This binary model will allow us to simplify the analytical solutions, although we may have significant errors when l ∼ s.

Fracture number density p 30
The p 30 measure counts the number of fractures N per unit volume V . Fracture occurrence can be measured as µ (l, s) = (l, s). The corresponding lacunarity λ p 30 under our hypothesis is then given by: If we consider that fractures have all the same size l f , then: If we consider now that the network follows a power-law size distribution with minimum and maximum fracture size l min and l max , the p 30 lacunarity equation is defined by three regimes: (1/p 30 ) · s −3 , s l max (12)

Fracture intensity p 32
Fracture intensity p 32 measures the total fracture surface per unit volume. The contribution of each fracture of size l in the domain S is µ (l, s) = (l, s) β f l 2 . The corresponding lacunarity λ p 32 is then given by: Following the same reasoning as for p 30 density, if we consider that all fractures have the same size l f , we find that the p 32 lacunarity λ p 32 follows Eq. (11). If we consider now that the network follows a power-law size distribution, the p 32 lacunarity equation is defined by three regimes:

Percolation parameter p
The percolation parameter gives an idea of the connectivity of the network (Bour andDavy, 1997, 1998;Robinson, 1983). Fundamentally, it is total excluded volume around fractures per unit volume so that for disk-shaped fractures the associated measure µ(l, s) = (π 2 /8) · (l, s) l 3 (De Dreuzy et al., 2000). The corresponding lacunarity λ p is then given by: If we consider that all fractures have the same size l f , then: If fracture size distribution follows a power-law size distribution, the percolation parameter lacunarity becomes:  . Theoretical (dash) and experimental (dots) p 10 lacunarity curves for fracture network sets (3-6) following power-law size distribution at bulk density p 32 = 0.5.

Density variability of DFN with power-law size distributions and various orientation distribution
In this section, we aim at validating the analytical solutions developed in Sect. 2, with numerical experiments on Discrete Fracture Networks (DFN). We generate some very simple Poissonian DFN models where the position of each fracture is set randomly within a cubic system of size L = 400, and other fracture parameters (orientations, size. . . ) are picked up in the corresponding distribution independently of each other. The fracture shapes are disks, and the diameters l are either constant or power-law distributed. In order to minimize finite-size effects, fracture centers are generated in a cubic box system of size (L + l max ) with l max the largest fracture. We here define two fracture network sets with constant size: (1) l f = 20 L and (2) l f = L = 400. We also define four fracture network sets (3, 4, 5, 6) following powerlaw size distribution with l min = 1, l max = 100, and where the power-law exponent a is within the range [2, 5]. All network sets are simulated for different fracture intensities p 32 ∈ {0.1, 0.2, 0.3, 0.4, 0.5}. We test two different fracture orientation distributions: (a) uniform, all orientations are equally represented, (b) Fisher distribution f (θ, κ) (Fisher, 1953) with a mean pole θ defined by a strike and dip angle both of 45 • , and a dispersion factor κ = 15. For statistical analysis, we perform 50 realizations of each fracture sets. Figure 1 shows examples of generated networks and their associated size distributions.

1-D measurements
We compute the p 10 lacunarity curves on 81 boreholes crossing the whole domain, equally spaced and parallel to the z direction, divided in subsamples of size s ∈ [0.02L, L] with no overlap. The general p 10 lacunarity curve of the system is obtained by stacking the contribution of each borehole. Figure 2 shows that the lacunarity curves associated to the generated Poissonian DFN models with power-law size distribution and mean density p 32 follow Eq. (9) as predicted by Darcel et al. (2013) and Lu et al. (2017). When s ∼ L, the lacunarity curve drops down because of finite-size effects. Lacunarity curves are not the same for the uniform and the Fisher orientation distributions, because the p 10 densities are not the same. Indeed, for the same p 32 density, the p 10 value measured on a borehole depends of the relative orientation between fractures and the borehole. For uniform orientation p 10 = p 32 /2, and for the defined Fisher orientation p 10 = 0.75p 32 which can be found from Wang (2005) formulae.

3-D measurements
To construct the experimental three-dimensional density lacunarity curves, we divide the cubic domain of size L in L 3 /s 3 subdomains with no overlap at different scales s. On each subdomain, the three-dimensional fracture densities defined in Sect. 2.2 (p 30 , p 32 , p) are computed numerically in order to obtain the associated mean and standard deviation over the whole domain at scale s, giving access to the corresponding lacunarities. For the analytical analysis, we consider that the fracture disk shape factor is β f = π/4, and that the subdomain shape factor is β s = 1. Figure 3 focuses on the p 32 lacunarity of constant size fracture networks with total fracture intensity p 32 = 0.5. For the case where the fracture size l f = 20, Eq. (10) shows two asymptotes for the experimental lacunarity curves for s l f and for s l f , respectively. When s ∼ l f , the equations overpredict the lacunarity by a factor of at most 3. Figure 4 focuses on various 3-D density lacunarity curves for networks following power-law size distributions. When the study scale s is larger than the minimum fracture size l min , the p 30 lacunarity curves (Fig. 4a) evolves as a threedimensional binomial process (∼ s −3 ), which is coherent with the fact that fractures are positioned randomly in space. The fracture intensity lacunarity λ p 32 and percolation parameter lacunarity λ p (Fig. 4b and c respectively), are divided in three different regimes. When we investigate the dependence upon the scale s ∈ [l min , l max ], their scaling clearly depends of the power-law size distribution exponent a. For large a values (a ≥ 5), the p 32 lacunarity curve scales as ∼ s −3 be-cause small fractures are dominant. For this range of study scale, we can notice a slight discrepancy between the analytical solutions and the numerical experiments. The lower a the larger this discrepancy. This can be explained by the fact that for low a values, the number of fractures whose size l ∼ s is always significant. Nevertheless, these analytical solutions reproduce well the observed scaling of fracture densities lacunarity. At scale s ∼ l max , the difference of p 32 lacunarity between networks following either a power-law size distribution of exponent a = 2 or a = 5 vary up to a factor 10 3 . Finally, Fig. 4d analyses λ p 32 at fixed scale s = 20 for different bulk p 32 values, highlighting that fracture density lacunarity is inversely proportional to the total fracture density. Finally, for any fracture density, all lacunarity curves are identical whatever the fracture orientation distribution (uniform or Fisher), which proves the orientation independency.

Conclusion
We propose here analytical solutions to quantify the density variability associated to Poissonian fracture networks, using the dimensionless variance parameter λ that we call lacunarity. Solutions are demonstrated for any kind of fracture density of any dimension, with application to uniform and power-law size distribution. The application to other fracture size distribution (exponential, log-normal. . . ) is straightforward. These analytical solutions were validated numerically by computing the density variability on Poissonian DFN models for different bulk density values, fracture size and orientation distributions. We show that the variability of three-dimensional densities is dependent on the study scale and the fracture size distribution, but not on the orientation distribution. Moreover, we show that for networks following power-law size distributions, the scaling of the variability of the fracture intensity p 32 is strongly dependent of this powerlaw exponent. This suggest that the variability of p 32 density cannot be estimated from borehole p 10 measurements, as it is done to quantify the mean p 32 density over the whole network using simple stereological rules. These solutions, which are developed for purely stochastic networks, can thus be used as a reference to estimate three-dimensional fracture density variability on the field, which is out of reach of our investigation technics. Further work is ongoing to quantify the fracture density variability for networks showing fractal correlations (Darcel et al., 2003b), which may also be quantified from the field. We expect the clustering effect associated to such networks to increase significantly the variability of fracture densities and its scaling.
Code availability. The Python script used for generation and analysis of Discrete Fracture Networks in this paper can be found at https://github.com/elavoine/DFNDensityVariability (last access: 30 May 2019).
Author contributions. EL, PD and CD conceived the idea of this study. EL and PD developed the analytical solutions. EL and RLG developed the code for generation and analysis of DFNs. EL took the lead in writing manuscript, all authors providing critical feedback.
Competing interests. The authors declare that they have no conflict of interest.
Special issue statement. This article is part of the special issue "European Geosciences Union General Assembly 2019, EGU Division Energy, Resources & Environment (ERE)". It is a result of the EGU General Assembly 2019, Vienna, Austria, 7-12 April 2019.