Jean-François Barthélémy
email: j-francois.barthelemy@ifp.fr

Martin L E Guiton

Jean-Marc Daniel

Estimates of fracture density and uncertainties from well data

Keywords: fracture networks, density probability law, scanline, cylindrical well, partial intersection

This paper aims at building a method to estimate the probability law governing the 3D fracture density of a fractured rock conditioned to the number of traces observed on a borehole image when the spatial distribution of fracture centers is assumed to follow a Poisson process. A closed-form expression of this law, allowing to calculate its mean value as well as a confidence interval, is derived in both cases of a lineic well (scanline) and a cylindrical well. The latter is better adapted to the situation of fracture size of the same order of magnitude as the well radius, which enables the presence of partial traces. In particular, the method takes into account the bias in the density estimate due to the fact that a fracture may cut the well along two distinct traces according to the considered fracture size. Monte Carlo simulations finally show a good agreement with the theoretical results of mean density and confidence interval.

Introduction

Naturally fracture networks in subsurface reservoirs have long been recognized to potentially influence the flow of hydrocarbon ( [START_REF] Stearns | Reservoirs in fractured rock. Stratigraphic oil and gas fields[END_REF], [START_REF] Nelson | Geological analysis of naturally fractured reservoirs[END_REF]), water or CO 2 ( [START_REF] Beliveau | Analysis of the waterflood response of a naturally fractured reservoir[END_REF], [START_REF] Beliveau | Analysis of a tertiary CO 2 flood in a naturally fractured reservoir[END_REF], [START_REF] Faure | Numerical simulation of CO 2 injection in a heavy oil fractured reservoir -Case study of Bati-Raman field, Turkey[END_REF]). At the scale of rock matrix heterogeneity, micro-fractures, which are the main ones fully or partially observed on core data, can provide a substantial contribution to the initial volume of fluid in place ( [START_REF] Narr | Estimating average fracture spacing in subsurface rock[END_REF], [START_REF] Howard | Description of natural fracture systems for quantitative use in petroleum geology[END_REF]). At a larger scale of a drainage area around a given well, natural fracture network may significantly enhance the permeability, as identified in pressure transient analysis with change of derivative in time, or production events like mud-loss or early water breakthrough. At the late stage of production, the presence of natural fractures often controls the efficiency of enhanced recovery processes ( [START_REF] Boerrigter | Fractured reservoir simulation: Case studies[END_REF], [START_REF] Faure | Numerical simulation of CO 2 injection in a heavy oil fractured reservoir -Case study of Bati-Raman field, Turkey[END_REF]). At shallow depth, natural fractures are also known to control the mechanical behavior of rock masses involved in the stability of civil engineering structures such as tunnels ( [START_REF] Barton | Deformation phenomena in jointed rock[END_REF], [START_REF] Goodman | Block theory and its application[END_REF]). Considering fractures as micro-heterogeneities in a representative volume element (R.V.E.) of the reservoir, their hydraulical or mechanical effects are involved in simulations often through homogenized properties. For instance, the macroscopic permeability as well as the percolation threshold can be estimated by numerical ( [START_REF] Koudina | Permeability of three-dimensional fracture networks[END_REF], [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF], [START_REF] Sarda | Hydraulic characterization of fractured reservoirs: simulation on discrete fracture models[END_REF], [START_REF] Lange | Hydraulic characterization of faults and fractures using a dual medium discrete fracture network simulator[END_REF], [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF], [START_REF] Garcia | Fast and efficient modeling and conditioning of naturally fractured reservoir models using static and dynamic data[END_REF]) or theoretical methods ( [START_REF] Berkowitz | Analysis of fracture network connectivity using percolation theory[END_REF], [START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 1-Effective connectivity[END_REF], [START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 2-Permeability[END_REF], [START_REF] Fokker | General anisotropic effective medium theory for the effective permeability of heterogeneous reservoirs[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF], [START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF]). Whatever the method, the construction of homogeneized property should require the following stages. First, it is necessary to collect information on parameters describing the fractures in the R.V.E. from scarce available data, namely mainly from core data analysis or from image log (e.g. FMS, FMI and UBI) interpretation ( [START_REF] Olarewaju | Stochastic simulation of fracture density for permeability field estimation[END_REF], [START_REF] Wu | Imaging 3-D fracture networks around boreholes[END_REF], [START_REF] Özkaya | Fracture length estimation from borehole image logs[END_REF]). Besides the qualitative description of a possible cementation or material infill which should be considered for conductivity, a quantitative description is generally first elaborated from geometrical parameters : orientation, size, spatial distribution, connectivity, opening and density. In the following, we will focus only on the fracture density (see definition in 1) which is particularly important to appropriately estimate homogenized hydraulical properties. The usual procedure is then to extrapolate fracture density from well data to the full field reservoir by means of a geostatistical approach, either with basic variogram analyses [START_REF] Olarewaju | Stochastic simulation of fracture density for permeability field estimation[END_REF] or with the construction of an explicative function relating fracture density to drivers available at the full field scale. For instance, multilinear explicative functions are obtained with discriminant analysis in [START_REF] Gauthier | Integrated fractured reservoir characterization: a case study in a North Africa field[END_REF], [START_REF] Gauthier | The characterization of a fractured reservoir: a multidisciplinary with application to an offshore, Abu Dhabi carbonate reservoir[END_REF] and [START_REF] Garcia | Fast and efficient modeling and conditioning of naturally fractured reservoir models using static and dynamic data[END_REF], while non-linear explicative functions are provided by neural network techniques in [START_REF] Ouenes | Practical application of fuzzy logic and neural networks to fractured reservoir characterization[END_REF] and [START_REF] Tran | Characterisation and modelling of naturally fractured reservoirs[END_REF]. The drivers used to explain fracture density may be interpreted according to the geological history during which natural fractures have formed. It is worth considering the influence of mechanical properties coming from lithology (e.g. porosity, grain size, dolomitization in [START_REF] Stearns | Reservoirs in fractured rock. Stratigraphic oil and gas fields[END_REF], [START_REF] Gauthier | Integrated fractured reservoir characterization: a case study in a North Africa field[END_REF]) or from sedimentology like the small bed thickness positively correlated to the spacing between joints ( [START_REF] Ladeira | Relationship between fracture spacing and bed thickness[END_REF], [START_REF] Narr | A method for estimating subsurface fracture density in core[END_REF], [START_REF] Narr | Joint spacing in sedimentary rocks[END_REF]). Natural fractures distribution may also have been controlled by tectonic deformation like the extension in fold outer-arc which correlates to strata curvature ( [START_REF] Lisle | Detection of zones of abnormal strains in structures using gaussian curvature analysis[END_REF], [START_REF] Hennings | Combining outcrop data and three-dimensional structural models to characterize fractured reservoirs: an example from Wyoming[END_REF], [START_REF] Gauthier | The characterization of a fractured reservoir: a multidisciplinary with application to an offshore, Abu Dhabi carbonate reservoir[END_REF]) and the large shear deformation at the vicinity of faults ( [START_REF] Pohn | Joint spacing as a method of locating faults[END_REF], [START_REF] Gauthier | Integrated fractured reservoir characterization: a case study in a North Africa field[END_REF]) which is often partly accommodated by pervasive fractures. A crucial point of the workflow is the link between the estimate of fracture density from well observations and the 3D density of a R.V.E.. Indeed, it is for instance well known that the first one is strongly influenced by the relative orientation of the well with respect to fractures, while R.V.E. density should not. Since [START_REF] Terzaghi | Sources of error in joint surveys[END_REF] and then in [START_REF] Dershowitz | Rock Joint System[END_REF], [START_REF] Dershowitz | Interpretation of fracture spacing and intensity[END_REF], [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF], [START_REF] Wang | Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces[END_REF] a correction factor is widely used to compensate the poor chance to sample fractures sub-parallel to the well trajectory. More recently, this method has been extended to account for the intersection of fractures of finite size with a well cylinder ( [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], [START_REF] Wang | Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces[END_REF]) or for the bed thickness ( [START_REF] Narr | A method for estimating subsurface fracture density in core[END_REF], [START_REF] Narr | Estimating average fracture spacing in subsurface rock[END_REF]). Alternatively, in [START_REF] Wu | Imaging 3-D fracture networks around boreholes[END_REF] fracture density is computed with data from different locations on the well but converging towards a same volume of density measurement. When a relationship has been derived between the mean density computed from the observations and that of a corresponding given 3D fracture network, the uncertainty of the density estimate remains to be quantified, in order to allow for instance the comparison of density estimate from different wells. An empirical way to quantify this uncertainty would be to realize many 3D discrete fracture networks (D.F.N.) in a Monte Carlo simulation with a priori postulate on the 3D model properties (orientation, size and density of fractures) following for instance the idea presented in [START_REF] La Pointe | Stochastic Estimation of Fracture Size Through Simulated Sampling[END_REF] in which the size distributions is constrained from length of observed fracture traces. Beside the obvious time consuming drawback of such a method, the quantification of a probability distribution for the 3D fracture density would lack of precision. Also, like in [START_REF] Dershowitz | Rock Joint System[END_REF], one could consider a distribution for the number of fractures in a given 3D volume, e.g. the discrete Poisson distribution (see details in 1), but to the authors' knowledge, the determination of the density parameter which characterizes such distribution given the well data remains to be done. In this context, this paper aims at proposing closed-form solutions to this issue. The method presented in this paper therefore allows to quantify efficiently the uncertainty of the density estimate and to envision a quality control on the sample data. After recalling the different definitions of the density and the hypothesis of Poisson distribution, the paper falls into two parts. The first one considers the well as a scanline which can be crossed by plane fractures and the second one explicitely takes into account the cylindrical geometry of the well on which partial fracture traces can be observed. For both models, we derive the complete probability law of the fracture density conditioned to the number of observed traces along a well segment of given length. Eventually, this law allows to calculate the mean fracture density as well as quantiles, namely the confidence interval for any given tolerance. Whereas no assumption is made on the fracture shape in the lineic case, the procedure proposed for the cylindrical well requires the knowledge of the size distribution of fractures assumed to be elliptic. This distribution has therefore to be estimated beforehand, for example by assuming a correlation with another fracture property which can be observed on borehole images, e.g. spacing or aperture [START_REF] Özkaya | Fracture length estimation from borehole image logs[END_REF], or by means of an inversion procedure constrained by the observation of fracture traces on outcrop surfaces (e.g. see [START_REF] Priest | Determination of discontinuity size distributions from scanline data[END_REF], [START_REF] Berkowitz | Stereological analysis of fracture network structure in geological formations[END_REF], [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF] and [START_REF] Tonon | Closed-form and numerical solutions for the probability distribution function of fracture diameters[END_REF] for the intersection between circular fractures and a plane outcrop, and see [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF] for the intersection between a 3D fracture network and a well cylinder in partial or full traces). However, according to the considered size distribution, the density probability presented in this paper takes into account the bias induced by the fact that two distinct traces can possibly come from the same fracture, which implies a reduction of the apparent density.

Hypotheses and definitions

As often assumed, the fractures will be considered as planar (see [START_REF] Zhang | Stereological relationship between trace length and size distribution of elliptical discontinuities[END_REF] for a review about this hypothesis). The fractured medium that is studied here may contain several sets of fractures. A set gathers fractures according to the distribution of normal vectors around a mean pole and a rather small dispersion around the pole. For example, a method to create fracture sets statisfying bivariate normal distributions can be found in [START_REF] Marcotte | Automatic joint set clustering using a mixture of bivariate normal distributions[END_REF] (other distributions are presented in [START_REF] Fisher | Statistical analysis of spherical data[END_REF]). Each set is also characterized by a size distribution. The presence of several sets may imply a correlation between the orientation and size distributions. Nevertheless, it is reasonable to assume these properties are independent within each set. In the following, a unique set is considered. It is then characterized by an orientation distribution, a size distribution as well as a density. The latter in a rock mass crossed by a well is quantified by means of the following definitions ( [START_REF] Dershowitz | Rock Joint System[END_REF], [START_REF] Dershowitz | Interpretation of fracture spacing and intensity[END_REF], [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF]): P 30 is the mean number of fracture centers per unit volume (in /m 3 ), P 32 is the mean fractured surface per unit volume (in m 2 /m 3 ), P 21 is the mean curvilinear length of fracture traces per well unit surface (in m/m 2 ), recalling that a trace is defined as the intersection between the fracture and the well surface (seen as a cylinder) and P 10 is the mean number of fracture traces per well unit length (in /m).

To identify proper parameters of the size, orientation and spatial distributions, it is necessary to consider studied domains, namely intervals, within which those parameters are stationary. Therefore, the intervals must be small enough to ensure the homogeneity of the distributions. The hypothesis of the independence between the size distribution and the spatial location is also argued in [START_REF] Zhang | Stereological relationship between trace length and size distribution of elliptical discontinuities[END_REF]. The last hypothesis, which will play a major role in the sequel, is that the fracture centers follow a Poisson distribution. Although this is not the only spatial distribution which has been identified on field observations, numerous works ( [START_REF] Priest | Discontinuity spacings in rock[END_REF], [START_REF] Rives | Joint spacing: analogue and numerical simulations[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF]) have proven its validity in a large number of field studies. In particular, it is shown in [START_REF] Rives | Joint spacing: analogue and numerical simulations[END_REF] that it may be well adapted to describe networks of low density or early stages of fracturation. The Poisson distribution means that the number N V of centers which can be found in a volume V is given by:

P NV (n V ) = e -P30 V (P 30 V ) nV n V ! = P P30V (n V ) (∀n V ≥ 0) (1)
where P λ (n) = e -λ λ n /n! denotes the discrete Poisson distribution of positive scalar parameter λ. It is very important to note that P 30 plays the role of a parameter of the law [START_REF] Andersson | A stochastic model of a fractured rock conditioned by measured information[END_REF]. The expectation of N V is equal to P 30 V , which gives P 30 its physical meaning as the mean number of centers per unit volume. Nevertheless, one must emphasize that P 30 , as an unknown parameter, should not be directly identified as the number of fracture centers divided by the volume of a given observation since the uncertainty on the actual value of P 30 would not be then taken into account, especially if the number of centers is too small. According to the previous definitions, this parameter P 30 can be replaced in (1) by the usual 3D density P 32 thanks to

P 32 = P 30 E [S] (2)
where E [S] is the expectation of the fracture surface S. The density P 32 now also appears as a positive scalar parameter allowing to write the Poisson distribution of the number of fracture centers in a given volume. One of the main issues of fracture network characterization is to identify P 32 . Many previous works ( [START_REF] Dershowitz | Rock Joint System[END_REF], [START_REF] Dershowitz | Interpretation of fracture spacing and intensity[END_REF], [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF], [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF], [START_REF] Wang | Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces[END_REF], [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]) consider P 32 as fixed albeit unknown. In addition, parts of those works are devoted to the derivation of relationships between different parameters: e.g. between P 32 and P 10 or between P 32 and P 21 . By this way, it is intended to identify P 32 by means of 1D or 2D samplings, for instance the number of fracture traces divided by the interval length in a scanline sampling. But this method may lose accuracy if the number of observations is small. By contrast, the present contribution aims at quantifying the uncertainty of the determination of the parameter P 32 from observations (scanline or cylindrical sampling).

2 Lineic data

Orientation and size distributions on a 1D well path

We consider here that the well trajectory is defined by a straight line. It is well known that any 3D statistics of a fracture property depending on the orientation and/or extension distributions of the fractures does not simply apply to a sampling of fractures intersecting the well path ( [START_REF] Terzaghi | Sources of error in joint surveys[END_REF], [START_REF] Dershowitz | Interpretation of fracture spacing and intensity[END_REF], [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF], [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF], [START_REF] Wang | Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces[END_REF]). Indeed, statistics are biased by the fact that the number of large fractures or the number of those which are the most perpendicular to the well path are over-estimated in such a sampling. Thus, the relationship between a 3D statistics and its 1D counterpart must involve a correction factor inversely proportional to the mean number of fractures of given orientation and surface crossing the well path. This number is proportional to the area of the projection of the fracture on the plane perpendicular to the well path. If Θ denotes the random acute angle between the fracture normal vector and the well direction (0 ≤ Θ ≤ π/2) and S the random surface of the fracture (see Fig. 1), the mean number of fractures such that θ ≤ Θ < θ + dθ and s ≤ S < s + ds intersecting a well segment of length L is:

dN L = P 30 L s cos θ p S (s) ds p Θ (θ) dθ (3)
where p S (resp. p Θ ) is the 3D p.d.f. (probability density function) of the variable S (resp. Θ). Consequently, the p.d.f. of the couple (S, Θ) defined on fractures intersecting the well path and denoted by p L S,Θ is then related to its 3D counterpart by:

p L S,Θ (s, θ) = s cos θ E [S cos Θ] p S,Θ (s, θ) = s E [S] p S (s) cos θ E [cos Θ] p Θ (θ) (4)
The second equality of (4) takes advantage of the independence of Θ and S 3D distributions. Moreover, (4) proves that Θ and S, as variables concerning only intersecting fractures, are also independent:

p L S,Θ (s, θ) = p L S (s) p L Θ (θ) with p L S (s) = s E [S] p S (s) and p L Θ (θ) = cos θ E [cos Θ] p Θ (θ) (5)
2.2 Density analysis: P 32 from P 10

Poisson process on a well path

The aim of this section is to prove that the fracture traces along the well path obey to a Poisson process if the 3D spatial distribution is also of Poisson type i.e. satisfies [START_REF] Andersson | A stochastic model of a fractured rock conditioned by measured information[END_REF]. For this purpose, we build a cylindrical volume around the well path as in [START_REF] Zhang | Estimating the intensity of rock discontinuities[END_REF] (see Fig. 1). The surface of the cylinder section is denoted by S c and the height by L. Assuming that S c is very large compared to the order of magnitude of the fracture surface, the probability for any fracture to cut the well path is:

L Θ S c
q = E [S cos Θ] S c = E [S] S c κ o with κ o = E [cos Θ] (6)
Let N ′ L be the number of fractures cutting the well path among those for which the centers belong to the cylindrical volume. It is worth observing that the probability for a fracture to have its center inside the cylindrical volume while cutting the well path is equal to the probability to have its center outside while cutting the well segment contained inside the cylinder. Hence, the probability law of the random number N L of fractures cutting this segment of length L is the same as that of N ′ L . Considering the definition of N ′ L , we deduce that its conditional probability given the number N V of centers inside the cylinder follows a binomial distribution:

P N ′ L |NV (n L |n V ) =      0 if n L > n V B q,nV (n L ) = C nL nV q nL (1 -q) nV -nL if n L ≤ n V (7)
where

C nL nV = n V !/(n L !(n V -n L )!)
is the binomial coefficient and q defined in [START_REF] Berkowitz | Analysis of fracture network connectivity using percolation theory[END_REF]. It is now possible to derive P NL by introducing [START_REF] Berkowitz | Stereological analysis of fracture network structure in geological formations[END_REF] and (1) in the marginal law formula [START_REF] Jourdain | Probabilitï¿ 1 2 s et Statistique[END_REF]

P NL (n L ) = P N ′ L (n L ) = +∞ nV =0 P N ′ L |NV (n L |n V ) P NV (n V ) = +∞ nV =nL B q,nV (n L ) P P30V (n V ) (8)
Exploiting ( 2), ( 6), V = S c L and the identity +∞ n=0

x n /n! = e x , (8) becomes:

PN L (nL) = e -P 32 κo L/q (P32 κo L) n L nL! +∞ n=0 P32 κo L (1 -q) /q n n! = e -P 32 κo L (P32 κo L) n L nL! = PP 32 κoL (nL) (9)
The probability law [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF] proves that the property of cutting the well path is a Poisson process. It follows that the mean number of traces per well unit length is consistent with the integration of (3) and highlights the correction factor between P 10 and P 32 already put in evidence in numerous works after [START_REF] Terzaghi | Sources of error in joint surveys[END_REF] ( [START_REF] Dershowitz | Interpretation of fracture spacing and intensity[END_REF], [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF], [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF], [START_REF] Wang | Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces[END_REF]):

P 10 = E [N L ] L = P 32 κ o (10)
Similarly to P 32 , P 10 plays the role of a parameter in the Poisson process characterizing the number of fracture traces on the well interval. Moreover the spacing ∆ between two consecutive traces obeys to an exponential law ( [START_REF] Dershowitz | Rock Joint System[END_REF]): p ∆ (δ) = P 10 e -P10 δ (∀δ ≥ 0)

As aforementioned, some field studies have put in evidence fracture spacings obeying to (11) ( [START_REF] Priest | Discontinuity spacings in rock[END_REF], [START_REF] Rives | Joint spacing: analogue and numerical simulations[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF]), which allows to define an observable criterion of the spatial distribution. Indeed, it is recalled in [START_REF] Odling | Variations in fracture system geometry and their implications for fluid flow in fractured hydrocarbon reservoirs[END_REF] that if the latter if of Poisson type, the ratio between the standard deviation and the mean of the spacing distribution is 1, whereas the ratio is smaller than 1 if the traces are more regularly spaced, and greater than 1 if the traces are clustered.

Probability law of P 32

The expression (9) has been obtained under the hypothesis that P 32 was a known parameter. Assuming now that P 32 is a random variable, (9) writes:

P NL|P32 (n L |p) = e -p κo L (p κ o L) nL n L ! (12)
Invoking Bayes' theorem [START_REF] Jourdain | Probabilitï¿ 1 2 s et Statistique[END_REF] by taking [START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 1-Effective connectivity[END_REF] as the conditional probability and considering the prior p P32 and the marginal P NL probabilities as normalizing constants, the p.d.f. of P 32 given the number N L is obtained as the posterior probability:

p P32|NL (p|n L ) = e -p κo L p nL (κ o L) nL+1 n L ! = Γ nL+1,κoL (p) (13)
where Γ nL+1,κoL denotes the p.d.f. of the gamma type [START_REF] Jourdain | Probabilitï¿ 1 2 s et Statistique[END_REF] depending on the parameters n L + 1 and κ o L.

In particular, the conditional expectation, mode and standard deviation can be deduced as:

E [P 32 |N L = n L ] = n L + 1 κ o L ; Mode [P 32 |N L = n L ] = n L κ o L ; Var [P 32 |N L = n L ] = √ n L + 1 κ o L (14)
It must be emphasized that the expectation and the mode ( 14) of P 32 given N L are not identical and, in particular, E [P 32 |N L = 0] = 0. The conditional p.d.f. ( 13) is represented on Fig. 2 for different values of L or N L for a given expectation ( 14)

E [P 32 |N L = n L ] = 1/m.
It quantitatively shows the influence of the length over which data are gathered and the number of data on the width of the p.d.f.. In particular, the ratio between the standard deviation and the expectation ( 14) decreases with respect to n L as 1/ √ n L + 1.

0 0.5 1.0 1.5 2.0 2.5 0 0.5 1.0 1.5 2.0 2.5 p p P32|NL (p|n l ) N L = 50, κ o L = 51 m N L = 40, κ o L = 41 m N L = 30, κ o L = 31 m N L = 20, κ o L = 21 m N L = 10, κ o L = 11 m Figure 2: P.d.f. p P32|NL for different values of L and N L such that E [P 32 |N L = n L ] = 1
The cumulative distribution function is:

F P32|NL (p|n L ) = P (P 32 ≤ p|N L = n L ) = 1 n L ! p κo L t=0 e -t t nL dt = 1 n L ! γ (n L + 1, p κ o L) ( 15
)
where γ is the so-called incomplete gamma function:

γ(a, x) = x t=0 e -t t a-1 dt (16)
Let us recall that the quantiles q α of the law of P 32 given N L are defined by:

F P32|NL (q α |n L ) = P (P 32 ≤ q α |N L = n L ) = α (17)
and can be obtained with respect to the integer n L by the resolution of:

q α = qα κ o L with 1 n L ! γ (n L + 1, qα ) = α (18)

Validation by Monte Carlo simulations

The validity of the P 32 p.d.f. found in the previous paragraph is now checked by means of Monte Carlo simulations. Discrete fracture networks are randomly drawn for several values of P 32 . The number of fractures N L intersecting a well segment, which length L is set to 20 in the following simulations, is calculated for each network. The cloud of points (N L /L, P 32 ) can then be compared on Fig. 3 to the expectation of the law [START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 2-Permeability[END_REF] as well as to the first and ninth deciles q 10 and q 90 respectively. For sake of simplicity, fractures are sometimes modeled as disks or more generally as ellipses characterized by two radii which follow scalar distributions such as power law, lognormal or exponential laws [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]. In the simulations of Fig. 3, the fractures are ellipses of aspect ratio 1/2 which largest radius follows a lognormal distribution of mean 1 m and standard deviation 0.3 m. As regards the orientation, a bivariate normal distribution [START_REF] Marcotte | Automatic joint set clustering using a mixture of bivariate normal distributions[END_REF] is adopted. It is recalled that such a distribution is characterized by an orthonormal frame (r, t 1 , t 2 ) (where r is the mean pole and (t 1 , t 2 ) is a basis of the plane tangent to the unit sphere at point r) or equivalently by three Euler angles (see Fig. 4) as well as by two standard deviation parameters s 1 and s 2 such that the orientation p.d.f. writes with respect to the unit sphere random vector n:

P (c i ≤ n • t i ≤ c i + dc i ; i = 1, 2) = 1 2 π s 1 s 2 I (s 1 , s 2 ) e -1 2 c 1 s 1 2 + c 2 s 2 2 dc 1 dc 2 (19)
where I (s 1 , s 2 ) is a normalization constant such that the integral of ( 19) over the unit disk c 2 1 + c 2 2 ≤ 1 is 1. For the simulations of Fig. 3, a mean pole inclined by θ = 45 1 2 with respect to the well direction e 3 has been chosen and s 1 = s 2 = 0.1. Moreover, the angle Ψ defining the direction of the major axis of a fracture in its plane (see Fig. 4) is taken constant and equal to π/2, i.e. the major axis belongs to the horizontal plane. As shown on Fig. 3, the law ( 13) is consistent with the Monte Carlo simulations; in particular, the proportion of points lying between the two decile curves among the 2000 drawn fracture networks is 79.16%.

3 Cylindrical data

Orientation and size distributions on a cylindrical well

In this paragraph, the well is considered as a circular cylinder of radius R c and the vector e 3 gives the direction of the axis. Because of the complexity of the intersection between a cylinder and any bounded surface representing a fracture, we limit ourselves to elliptic fractures of random major radius A and small radius B. As discussed in [START_REF] Zhang | Stereological relationship between trace length and size distribution of elliptical discontinuities[END_REF], the elliptical shape, albeit simple, is of great interest because of its ability to account for anisotropy in the fracture plane unlike the classically assumed circular shape To achieve the determination of such a fracture in the 3D space, we also have to know the center coordinates

(x o 1 , x o 2 , x o
3 ) and the angles (Θ, Φ, Ψ) defining the local frame of the ellipse: Θ still denotes the angle between the normal vector and the well direction, Φ is the angle between the projection of the normal vector onto the plane perpendicular to the well and the x 1 axis (directly related to the dip-azimuth if the well path is vertical) and Ψ allows to determine the rotation of the fracture in its own plane (see Fig. 4). The projection of such an ellipse onto the plane orthogonal to the well direction is an ellipse defined thanks to the quadratic form:

   x 1 -x o 1 x 2 -x o 2    • R • H -1 • R -1 •    x 1 -x o 1 x 2 -x o 2    ≤ 1 (20)
R =    cos Φ -sin Φ sin Φ cos Φ    and H =    cos 2 Θ A 2 cos 2 Ψ + B 2 sin 2 Ψ A 2 -B 2 cos Θ cos Ψ sin Ψ A 2 -B 2 cos Θ cos Ψ sin Ψ A 2 sin 2 Ψ + B 2 cos 2 Ψ    ( 21
)
The radii of the ellipse  and B ≤  in the (x 1 , x 2 ) plane are the square roots of the eigenvalues of H (21) and thus do not depend on Φ:

 (A, B, Θ, Ψ) = tr H + (tr H) 2 -4 det H 2 ; B (A, B, Θ, Ψ) = tr H -(tr H) 2 -4 det H 2 (22)
The product of the two radii, which could also have been calculated from the projection of the surface, is:

 B = A B cos Θ (23)
In the case of an initial circular fracture namely A = B, the simple expressions of  and B are obviously found:

 = A ; B = A cos Θ (24)
As shown by [START_REF] Gauthier | The characterization of a fractured reservoir: a multidisciplinary with application to an offshore, Abu Dhabi carbonate reservoir[END_REF], the orientation of the ellipse is given by a rotation of angle Φ of the eigenvectors of H. The general expressions [START_REF] Goodman | Block theory and its application[END_REF] and the determination of the axis of the projected ellipse show that the shape and the orientation of the latter are not independent whereas it is the case for the ellipse in the 3D space.

Statistics on fractures cutting the well

In the sequel, we follow the reasoning presented in [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF] to calculate the surface of the domain to which the projection of the fracture center must belong so that the fracture cuts the cylinder, given the radii (A, B) and the angles (Θ, Ψ). As for the 1D case, this surface, denoted hereafter by σ, will be very useful to define the p.d.f. of any variable over fractures intersecting at least partially the cylinder surface. The mean number of fractures such that θ ≤ Θ < θ + dθ, ψ ≤ Ψ < ψ + dψ, a ≤ A < a + da, and b ≤ B < b + db, intersecting a part of the well cylinder of length L is then:

dN L = P 30 L σ(a, b, θ, ψ) p A,B (a, b) da db p Θ,Ψ (θ, ψ) dθ dψ (25)
From the same reasoning as that leading to (4) in 1D, we write the p.d.f. of (A, B, Θ, Ψ) defined on fractures intersecting the well cylinder denoted by p C A,B,Θ,Ψ from the 3D p.d.f. p A,B and p Θ,Ψ :

p C A,B,Θ,Ψ (a, b, θ, ψ) = σ E [σ] p A,B (a, b) p Θ,Ψ (θ, ψ) (26)
where

E [σ] = a,b,θ,ψ σ(a, b, θ, ψ) p A,B (a, b) p Θ,Ψ (θ, ψ) dθdψdadb (27)
The expression of σ is now presented from results in [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF] and [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]. Actually σ depends on (A, B, Θ, Ψ) through ( Â, B) given by [START_REF] Goodman | Block theory and its application[END_REF]. Let us introduce the angle Λ such that

B = Â cos Λ (0 ≤ Λ ≤ π/2) (28)
It is worth noting that if A = B then Λ = Θ because of [START_REF] Hennings | Combining outcrop data and three-dimensional structural models to characterize fractured reservoirs: an example from Wyoming[END_REF]. The expression of σ derived in [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF] and recalled in [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF] can be simplified using classical elliptic integrals (details can be found in appendix A):

σ =            π (cos Λ Â2 + R 2 c ) + 4 E (sin Λ)  R c if  ≥ R c 8 E (sin Λ)  R c if  ≤ R c cos Λ 8 E (sin Λ)  R c + J -(η) if R c cos Λ <  < R c (29)
where E(k) is the complete elliptic integral of the second kind recalled in (64) (it can be estimated by

E (sin Λ) ≈ π/2 (1 + cos 2 Λ)/2 [51] if Λ is not in the neighborhood of π/2
), J -is defined in (63) and η in (66). The exact expressions [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF] require to resort to tabulated values of elliptic integrals allowing then to perform any numerical integration scheme involving σ such as the trapeze method, gaussian quadrature etc., to compute for example [START_REF] Jourdain | Probabilitï¿ 1 2 s et Statistique[END_REF]. Nevertheless, some authors ( [START_REF] Özkaya | Fracture length estimation from borehole image logs[END_REF], [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]) propose simpler approximate expressions valid only for Λ far enough from π/2 and which can be useful in some particular cases of exact integration:

σ ≈ σ a =      π (  + R c ) ( B + R c ) if  > R c 2 π (  + B) R c if  ≤ R c (30)
Unlike the 1D case, it must be emphasized here that the variables of shape and orientation concerning fractures intersecting the well surface are not in general independent. Indeed, the p.d.f. ( 26) is not necessarily the product of each marginal p.d.f. since σ, depending on (A, B, Θ, Ψ) through ( Â, B), is not the product of a function of (A, B) on the one hand and a function of (Θ, Ψ) on the other hand in [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF] or [START_REF] La Pointe | Stochastic Estimation of Fracture Size Through Simulated Sampling[END_REF]. Nevertheless, let us cite two cases for which the couples (A, B) and (Θ, Ψ) remain independent. First, if A ≥ B ≫ R c and Θ is far enough from π/2 then, using ( 23) and ( 28), [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF] implies

σ ≈ π Â B = S cos Θ (31)
This means that the p.d.f. ( 26) is the same as that found in the 1D case (4), which can be explained by the fact that the well radius is so small compared to the fracture length that the well can be considered as a line [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]. The second case is obtained when A < R c cos Θ and the fractures are circular (i.e. Â = A = B and Λ = Θ). Under these assumptions, the surface σ given by ( 29) proves then that A and Θ are independent:

σ = 8 A R c E (sin Θ) (32)

Fractures with full intersections

Instead of considering the population of fractures partially intersecting the cylinder, it could be interesting to focus on those characterized by a full intersection (this set is not empty if the condition R c < B is fulfilled), inducing then the definition of a new domain of surface σ f to which must belong the projection of the fracture center. The expression of σ f presented in [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF] does not seem to be right. Indeed, from calculations detailed in appendix A, it comes that:

σ f =                    0 if  cos Λ ≤ R c else J -(π/2) = π (cos Λ Â2 + R 2 c ) -4 E (sin Λ)  R c if  cos 2 Λ ≥ R c J -(ξ) if  cos 2 Λ < R c (33)
where the function J -is defined by (63) and ξ by (65). From the same reasoning as that leading to [START_REF] La Pointe | Stochastic Estimation of Fracture Size Through Simulated Sampling[END_REF], when Λ is far from to π/2, this surface can also be roughly estimated by ([44], [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]):

σ f ≈ σ f a =      π ( Â -R c ) ( B -R c ) if R c < B ≤ Â 0 if B ≤ R c (34)
As mentioned in [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF], the probability for a fracture cutting the well to fully cut it, given its orientation and length, is simply the ratio σ f /σ. Eventually, given size and orientation distributions, the proportion of fractures fully cutting the cylindrical well among the entire set cutting the well i.e. the probability that any trace be a full trace is then:

P (trace is full) = E σ f E [σ] (35)
This result is of major interest because it can readily be compared to the empirical ratio observed on the borehole image (full trace number over total trace number). Recalling that the orientation distribution can be estimated from well traces in an easier way than the size one, [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF] can somehow be used to assess the likelihood of the latter taken as an input of the model. Nevertheless, numerical simulations in [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF] showed that, in some cases, two different fracture size distributions could lead to very close trace length distributions. To improve the determination of the fracture size law which is assumed to be known for our density analysis, some specific methods can be found in the literature. For example, considering fractures as disks, an inversion method from the trace histogram on a plane outcrop is provided in [START_REF] Berkowitz | Stereological analysis of fracture network structure in geological formations[END_REF] and can be used for a cylinder if the fracture radius is much lower than the well radius. Another method, presented in [START_REF] Özkaya | Fracture length estimation from borehole image logs[END_REF], is based on a linear correlation between the fracture size and an observable property (spacing, aperture etc.) fitted by comparison between [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF] and its empirical counterpart.

Fractures with single or double traces

A nonempty intersection between a cylinder and an ellipse can be either a single trace (partial or full) or a double trace. The latter case occurs when the boundary of the ellipse or the projection of the ellipse onto the plane perpendicular to the cylinder axis cuts four times the cylinder. The intersections can be practically found, for instance for simulations, by inserting in [START_REF] Gauthier | The characterization of a fractured reservoir: a multidisciplinary with application to an offshore, Abu Dhabi carbonate reservoir[END_REF] the following parameterization of the well circle

x 1 (t) = R c 1 -t 2 1 + t 2 ; x 2 (t) = R c 2 t 1 + t 2 (36)
The equation obtained by replacing the inequality (20) by an equality is a fourth-degree polynom in t. The number of real roots determines the number of traces : four reals roots imply two traces, two real roots imply one trace and no real root corresponds either to an empty intersection between the cylinder and the ellipse or to a full intersection depending on the sign of the left hand side of (20) applied on any point of the circle [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF]. The surface σ 2 of the domain to which the center of the projected ellipse must belong in order to show two traces on the well can be calculated by the procedure given in appendix A leading to:

σ 2 =                            0 if  cos 2 Λ ≥ R c or  ≤ R c cos Λ else -J -(π/2) if  cos Λ ≤ R c ≤  -J -(η) if  < R c J -(ξ) -J -(π/2) if  cos Λ > R c (37)
where J -is expressed in (63) (J -(π/2) in ( 33)), η in (66) and ξ in (65). Knowing σ (29) and σ 2 [START_REF] Narr | Estimating average fracture spacing in subsurface rock[END_REF], it is now straightforward to calculate the surface σ 1 of the domain in which the center of the projected ellipse must be located so as to generate a single trace on the well:

σ 1 = σ -σ 2 (38)
3.2 Density analysis: P 32 from P 10

As previously put in evidence, the main difference with the lineic case is that the correction factor involved when sampling on fractures intersecting the cylindrical well depends on the fracture size distribution and not only on the orientation distribution. Whereas the latter can be rather easily estimated by fitting with sinusoids the fracture traces observed on borehole images, this is not the case concerning the size distribution.

As previously mentioned, methods for estimating the latter can be found in the literature ( [START_REF] Berkowitz | Stereological analysis of fracture network structure in geological formations[END_REF], [START_REF] Özkaya | Fracture length estimation from borehole image logs[END_REF], [START_REF] Priest | Determination of discontinuity size distributions from scanline data[END_REF], [START_REF] Tonon | Closed-form and numerical solutions for the probability distribution function of fracture diameters[END_REF]). However this work does not aim at determining the size distribution but at estimating the fracture density probability, in particular the mean value and uncertainties, knowing the orientation and size distribution. The latter should come from hypotheses and from a preliminary work which is the subject of current research.

Poisson process on a cylindrical well

Since a fracture can leave two traces on the well, it is necessary to distinguish the number N L of fractures which cut the well (fully or partially with one or two traces) and the number N t L of observed traces. A fracture leaving a double trace on the well contributes only once in N L but twice in N t L . Thus, denoting the fracture numbers having a single and a double trace respectively by N 1 and N 2 , it comes:

N L = N 1 + N 2 ; N t L = N 1 + 2 N 2 (39)
As regards the cylindrical data, the notion of trace within a cylinder of given length (studied interval) has to be clarified since some traces can have a part outside the studied interval. We adopt the convention to relate a trace to the studied interval if the latter contains the mid-point of the trace, which is defined as the point of curvilinear abscissa (s A + s B )/2 where s A and s B are the curvilinear abscissae of the extremities, possibly outside the studied interval. Obviously this definition requires to know the whole trace. The aim of this paragraph is to find the probability of the observation namely the number of traces N t L with respect to the fracture density P 32 considered here as a parameter of the Poisson distribution. In the case of a very small well radius compared to the fracture length, this task, developped in section 2.2, is rather easy because the number of traces is equal to the number of fractures cutting the well and, as the proportion of full intersections asymptotically reaches 100%, the size distribution is not involved in the probability [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF]. In the present case, the task is made more complex since the relationship between the observation and the fracture density may be biased by the presence of double traces and influenced by the fracture size. E [N L ] can be found from the knowledge of the orientation and size distributions by integration of (25) and a similar reasoning leads to E [N i ] i=1,2 :

E [N L ] = P 30 L E [σ] ; E [N i ] = P 30 L E [σ i ] (i = 1, 2) (40)
Eventually, (2) allows to write the expectations [START_REF] Nelson | Geological analysis of naturally fractured reservoirs[END_REF] with respect to P 32 :

E [N L ] L = P 32 κ with κ = E [σ] E [S] ; E [N i ] L = P 32 κ i with κ i = E [σ i ] E [S] (i = 1, 2) (41)
The ratios κ and κ i play the same role as κ o = E [cos Θ] in [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF]. The reasoning of section 2.2.1 leading to the probability law ( 9) can be adapted to the present case, observing that the number of fractures which centers belong the control cylinder of Fig. 1 and cutting the well follows the same probability law as the number of fractures leaving a trace on the studied cylinder part of length L whatever the position of the center. Following then the reasoning of section 2.2.1, it comes out that the numbers of fractures intersecting the well or having a single trace or a double trace follow Poisson laws:

PN L (nL) = PP 32 κL (nL) = e -P 32 κ L (P32 κ L) n L nL! and PN i (ni) = PP 32 κ i L (ni) = e -P 32 κ i L (P32 κi L) n i ni! (42)
Obviously these results can also be transposed to the number of full intersections using

κ f = E [σ f ] /E [S]
in the Poisson law. Considering the different combinations of single and double trace intersections, the probability law of the trace number then writes:

P N t L n t L = ⌊n t L /2⌋ n 2 =0 PN 1 n t L -2 n2 PN 2 (n2) = ⌊n t L /2⌋ n 2 =0 PP 32 κ 1 L n t L -2 n2 PP 32 κ 2 L (n2) (43)
where ⌊x⌋ denotes the integer part of the positive real x. It is useful to remark that κ = κ 1 + κ 2 thanks to [START_REF] Narr | A method for estimating subsurface fracture density in core[END_REF] and ( 41) and to introduce the probabilities to have a double or a single trace for a fracture cutting the well respectively denoted by r and 1 -r:

P (2 traces) = E [σ 2 ] E [σ] = κ 2 κ = r ; P (1 trace) = E [σ 1 ] E [σ] = κ 1 κ = 1 -r (44)
Indeed, after simple algebraic manipulations, (43) can be rewritten

P N t L n t L = ⌊n t L /2⌋ n2=0 B r,n t L -n2 (n 2 ) P P32κL n t L -n 2 (45)
where B r,n t L -n2 denotes the binomial probability of parameters r and n t L -n 2 already introduced in (7).

Probability law of P 32

A first simple approach to analyze P 32 is to focus only on full intersections observed on borehole images. In this case, all the reasoning presented in section 2.2.2 is still valid provided that κ o be replaced by κ f and that the number of traces be understood as the number of full traces. Nevertheless, considering cases for which the fracture length are of the same order of magnitude as the well radius, full traces may not exist or at least be so few that the P 32 estimate would not be reliable. Consequently, it is worth providing a P 32 estimate based on partial intersection analysis.

The probability law of P 32 in the cylindrical case follows from (45) by analogy with the reasoning allowing to deduce (13) from ( 12) namely as a posterior probability obtained by normalizing [START_REF] Pohn | Joint spacing as a method of locating faults[END_REF] thanks to the identity

⌊n t L /2⌋ n2=0 B r,n t L -n2 (n 2 ) = (1 -(-r) n t L +1 )/(1 + r) p P32|N t L p|n t L = 1 + r 1 -(-r) n t L +1 ⌊n t L /2⌋ n2=0 B r,n t L -n2 (n 2 ) Γ n t L -n2+1,κL (p) (46)
with Γ n,λ (p) = e -λp λ n p n-1 /(n -1)! already introduced in (13). The expectation and standard deviation are then given by:

E P 32 |N t L = n t L = n t L + 1 -f (n t L , r) κ L ; Var [P 32 |N t L = n t L ] = n t L + 1 -f (n t L , r) κ L (47)
with

f (n t L , r) = r n t L (1 + r) 1 + (-r) n t L -(1 -r) 1 -(-r) n t L (1 + r) 2 1 -(-r) n t L +1 (48)
The expressions [START_REF] Priest | Discontinuity spacings in rock[END_REF] differ from ( 14) by the term f (n t L , r) (48) which appears as a correction due to the presence of double traces. Indeed, in the absence of the latters (r = 0) or if they are simply just ignored, meaning that the traces are supposed to come from different fractures, this correction term vanishes. Ignoring double traces when r is far from 0 could lead to an important error in the P 32 estimate. Neglecting the terms r n t L for a large number of observations and r < 1, (48) and ( 47) become:

f (n t L , r) ≈ n t L →∞ r 1 + r n t L - r (1 -r) (1 + r) 2 (49)
and

E P 32 |N t L = n t L ≈ n t L →∞ 1 κ L n t L 1 + r + 1 + 3 r (1 + r) 2 ; Var P 32 |N t L = n t L ≈ n t L →∞ 1 κ L n t L 1 + r + 1 + 3 r (1 + r) 2 (50)
Keeping only the term in n t L for large values, the expressions [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF] show that the ratio between the standard deviation and the expectation still varies as 1/ n t L and that the expectation and standard deviation could be deduced from ( 14) by a correction on the trace number of 1/(1 + r) = (κ 1 + κ 2 )/(κ 1 + 2κ 2 ). This term exactly corresponds to the ratio between the mean fracture number intersecting the well and the mean trace number taking into account the double traces. This asymptotical behavior is obviously not valid for a small number of observed traces along L. The determination of quantiles as defined in ( 17) requires now to derive the cumulative distribution function associated with [START_REF] Priest | Determination of discontinuity size distributions from scanline data[END_REF]:

F P32|N t L p|n t L = 1 + r 1 -(-r) n t L +1 ⌊n t L /2⌋ n2=0 B r,n t L -n2 (n 2 ) (n t L -n 2 )! γ n t L -n 2 + 1, p κ L (51)
with the incomplete gamma function introduced in ( 16). Unlike the lineic case, the quantiles can not be simply deduced from tabulations of γ. Nevertheless, it is still possible to normalize them

q α = qα κ L with 1 + r 1 -(-r) n t L +1 ⌊n t L /2⌋ n2=0 B r,n t L -n2 (n 2 ) (n t L -n 2 )! γ n t L -n 2 + 1, qα = α ( 52
)
so that qα = qα (r, n t L ) do not depend on κ and L. The normalized quantiles can be searched by means of an iterative method such as a Newton-Raphson scheme:

q0 α = E P 32 |N t L = n t L ; κ L = 1 and qi+1 α = qi α - F P32|N t L qi α |n t L ; κ L = 1 -α p P32|N t L (q i α |n t L ; κ L = 1) (53)
Let us now consider a study case in which the well radius is R c = 0.1 m, the analyzed length is L = 20 m, and the number of observed traces is N t L = 20. The major radius A of the ellipse follows a lognormal distribution of mean 1 m and standard deviation 0.1 m and the aspect ratio is 0.3. The distribution of the fracture normal is of bivariate normal type [START_REF] Garcia | Fast and efficient modeling and conditioning of naturally fractured reservoir models using static and dynamic data[END_REF] such that the mean pole is inclined by θ = 80 1 2 with respect to the well direction (the fractures are almost parallel to the well axis) and the deviations are s 1 = s 2 = 0.1. The major axis of the ellipse belongs to the plane orthogonal to the well (Ψ = π/2 on Fig. 4). After computations of E [σ] (27) and E [σ i ] performed by means of gaussian quadrature with ten points along each dimension of integration, the previous data lead to a ratio r = 34.4% and to κ = 0.174 allowing then to draw the p.d.f. of Fig. 5. The discrepancy between the two curves allows to quantify the error made by neglecting the double traces in the reasoning, in other words by considering that the number of traces is equal to the number of fractures intersecting the cylinder. The normalized histogram also presented on Fig. 5 corresponds to numbers of discrete fracture networks generated in the validation process of section 3.2.3 and satisfying the same number of traces (here 20) for different P 32 intervals. This empirical p.d.f. shows a good agreement with the theoretical p.d.f. ( 46) taking into account the double trace correction.

Validation by Monte Carlo simulations

Similarly to the lineic case in section 2.2.3, the conditional p.d.f. of P 32 given a trace number on a segment of the cylindrical well is compared to Monte Carlo simulations. In particular, the technique allowing to take into account the effect of double traces on the density has to be confronted with simulations. That is why the size and orientation numerical parameters of the considered fracture set are chosen identical to those used in section 3.2.2 for which the mean proportion of fractures having a double trace is 34.4%. About 2000 fracture networks are drawn for different values of P 32 such that the maximum number of traces is 24 along a segment of length L = 20 i.e. 0 ≤ N t L /L ≤ 1.2. For a given fracture network, the number of (full, partial or double) traces is computed by means of the procedure presented in section 3.1.3. The corresponding couples (N t L /L, P 32 ) are plotted on Fig. 6 as well as the expectations and the 1 st and 9 th deciles of P 32 given the number of traces per unit length N t L /L according to three models: (M 1) this model takes into account the effect of double traces, hence the expectation is calculated by [START_REF] Priest | Discontinuity spacings in rock[END_REF] and the quantiles by ( 52), (M 2) this model ignores the possibility of having double traces and thus considers a bijection between the set of traces and the set of intersecting fractures (this amounts to take r = 0 in ( 47) and ( 52)) but still takes into account the effect of the fracture size distribution through the coefficient κ, (M 3) this model is based on the lineic results ( 14) and [START_REF] Fokker | General anisotropic effective medium theory for the effective permeability of heterogeneous reservoirs[END_REF] and differs from M 2 by the fact that the correction due to the fracture size is neglected.

The proportion of points lying between the two decile curves of the M 1 model is 80, 96% and only 60% between the curves of M 2. Moreover, Fig. 6 shows again the error made by neglecting the decrease of density estimate due to the presence of double traces. Nevertheless it is worth recalling that, besides the assumption of an elliptic shape of the fractures, M 1 and M 2 models require the knowledge of the fracture size distribution. Although the issue of relating the trace length distribution to that of the fractures has already been solved in [START_REF] Berkowitz | Stereological analysis of fracture network structure in geological formations[END_REF] in the case of circular fractures intersecting a plane, it may not be so obvious to transpose this result to a cylindrical well. A numerical treatment of this issue is the subject of a current work. But, without any information about the fracture size, the only available model is M 3, which only takes into account the correction due to the observed orientations as if the traces were full, generated by very large fractures. Obviously, the Fig. 6 shows that this model can not be satisfactory for the fracture set considered here. Indeed, the correction due to the fact that the fractures are almost aligned with the well leads to considerably overestimate the actual P 32 density. 4 Application: uncertainty on the fracture number in a D.F.N.

As already stated, the knowledge of the density, among other characteristics of the fracture network, allows to study the percolation and transport properties of fractured rocks either by numerical ( [START_REF] Koudina | Permeability of three-dimensional fracture networks[END_REF], [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF], [START_REF] Sarda | Hydraulic characterization of fractured reservoirs: simulation on discrete fracture models[END_REF], [START_REF] Lange | Hydraulic characterization of faults and fractures using a dual medium discrete fracture network simulator[END_REF], [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF], [START_REF] Garcia | Fast and efficient modeling and conditioning of naturally fractured reservoir models using static and dynamic data[END_REF]) or theoretical methods ( [START_REF] Berkowitz | Analysis of fracture network connectivity using percolation theory[END_REF], [START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 1-Effective connectivity[END_REF], [START_REF] De Dreuzy | Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 2-Permeability[END_REF], [START_REF] Fokker | General anisotropic effective medium theory for the effective permeability of heterogeneous reservoirs[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF], [START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF]). Thus, the uncertainty on the density should have repercussions on any subsequent calculations depending on this density. To illustrate the consequence of the uncertainty on P 32 , let us consider the problem of the generation of a D.F.N. consistent with well data and more precisely the number of fractures N V in a given volume V . In order to create the D.F.N. corresponding to a given fracture set, it is necessary to give the size and orientation distributions, whatever how they have been calibrated or even postulated. As already recalled, well data can be used to identify the orientation distribution whereas they are often insufficient to calibrate the size and outcrop data may be needed. Provided that these two distributions be known, the 3D density is still needed and it is assumed here that it has to be deduced from a scanline sampling of n L traces observed over an interval of length L.

As shown in section 2.2.2, P 32 is unknown and thus must be considered as a random variable following the law (13) depending on the correction factor κ o (10) introduced in (6). The probability law of N V given N L is then obtained by the following marginal law:

P NV |NL (n V |n L ) = +∞ p=0 P NV |P32 (n V |p) p P32|NL (p|n L ) dp (54)
It is worth noting that (54) rigorously applies if the observed fractures do not have to be taken into account in the number N V . Otherwise, P NV |P32 should be replaced by the probability of N V conditioned to P 32 as well as N L . This conditional probability would then differ from the Poisson distribution because N L fractures among N V are already known [START_REF] Andersson | A stochastic model of a fractured rock conditioned by measured information[END_REF]. Nevertheless, if the generation volume size is such that N V is large enough with respect to N L , these two laws are only slightly different. Recalling that the probability of N V given P 32 is of Poisson type (1) where P 30 is replaced by (2), the law of N V given N L writes after calculation:

P NV |NL (n V |n L ) = C nL nV +nL V nV (κ o L E [S]) nL+1 (V + κ o L E [S]) nV +nL+1 (55)
It appears on (55) that the conditional law of N V given the trace number on a scanline interval is not of Poisson type. The expectation and standard deviation of such a discrete law are:

E [N V |N L = n L ] = (n L + 1) V κ o L E [S] ; Var [N V |N L = n L ] = (n L + 1) V (V + κ o L E [S]) κ o L E [S] (56)
By contrast, the point of view consisting in neglecting the random character of P 32 and choosing an estimate p est as if it was deterministic obviously makes N V follow the initially assumed Poisson law for which the expectation and standard deviation are given by:

E N V |P 32 = p est = p est V E [S] ; Var [N V |P 32 = p est ] = p est V E [S] (57)
For instance, p est could be chosen as the empirical mean of P 10 corrected by κ o [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF] namely

p est = n L /(κ o L)
or it could also be estimated by the expectation ( 14) namely p est = (n L + 1)/(κ o L). The ratio between the standard deviation and the expectation in (57) varies in 1/ √ V , which simply recalls that the uncertainty on N V relatively to its mean value decreases when the generation volume grows in a Poisson process considering P 32 as a fixed value. Moreover, it is pointed out in [START_REF] Dershowitz | Rock Joint System[END_REF] that, for large density values, the ratio N V /V tends towards a continuous variable following a normal distribution. On the contrary, if the random character of P 32 is taken into account according to the fact that the only observation is N L = n L , the ratio between the standard deviation and the expectation given by (56) tends towards a finite value 1/ √ n L + 1 when V tends towards infinity. This means that, for large values of V , an increase of the latter does not significantly reduce the uncertainty on N V . As expected, the only way to reduce this uncertainty is then to get more observed data. By the way, the number of observed traces which is necessary to reach a given ratio ρ between the standard deviation and the expectation is:

n L = 1 + κ o L E [S] V ρ -2 -1 (58)
The two laws coming from the two approches to characterize the probability of N V are presented on Fig. ?? with numerical values which are representative of a common case of petroleum reservoir engineering (10 observed traces on 10 m, κ o = 1, E [S] = 10 m 2 and V = 4000 m 3 ): the first one (solid curve) conditioned to the observation N L and following the law ( 55) and the second one considering the two aforementioned fixed values of P 32 (dashed curves) and following a Poisson distribution. This figure clearly highlights the additional uncertainty introduced by the fact that P 32 is not precisely calibrated by a small number of traces on a well interval. Indeed, it appears that, to be consistent with the fixed well data, the range of possible values for N V is very large (between about 200 and 700 fractures) and should be taken into account for subsequent calculations as percolation thresholds or effective permeabilities. By contrast, the two other curves show a reduction of the range of possible values of N V by considering that P 32 is fixed and this relative reduction would even be more obvious for larger values of V . Hence, the model with fixed P 32 could prevent the user from retaining values of N V which could be consistent with the well data.

Conclusions and further works

In this paper, the probability law of the P 32 density conditioned to the number of traces observed within a well log interval has been derived. The model is based under the following assumptions: fractures are planar, the spatial distribution of their centers is of Poisson type and the orientation and size distributions are independent. Moreover, the studied well interval is assumed to be such that the spatial, size and orientation distributions can be considered as homogeneous. In the case of a lineic model of well, it is shown that the P 32 law is of the gamma type. The latter depends only on the number of traces, the length of the analyzed well segment and a correction factor which takes into account the probability law of the inclination of fractures with respect to the well direction. The second model considers the well as a cylinder on which partial fracture traces can be observed and assumes that fractures are elliptical. This model leads to a probability law of the P 32 density which writes as a linear combination of gamma probability density functions depending on the number of traces, the length of the analyzed well segment as well as correction factors. These factors depend not only on the fracture orientation but also on the size distribution unlike the lineic case. Indeed, they require the calculation of surfaces of domains to which the projection of a fracture center must belong to intersect the cylinder along one or two traces. Explicit relationships are provided between those surfaces and a priori assumed fracture size. It is believed that the present closed-form results should improve the efficiency of methods estimating the 3D fracture density from well data. In particular, one could readily compute the p.d.f. of the P 32 measure of density from a P 10 value of a moving window log, without the need to generate discrete fracture networks in Monte Carlo simulations. In complement to the results of many previous works related to this issue ( [START_REF] Dershowitz | Rock Joint System[END_REF], [START_REF] Dershowitz | Interpretation of fracture spacing and intensity[END_REF], [START_REF] Mauldon | Intersection probabilities of impersistent joints[END_REF], [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF], [START_REF] Wang | Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces[END_REF], [START_REF] Gupta | Stereological analysis of fracture networks along cylindrical galleries[END_REF]), the present work adds a contribution to the quantification of the estimate uncertainty by providing a probability density function of the 3D Poisson model density parameter given the well observation.

The proposed method requires the preliminary knowledge of the orientation distribution of the fractures which can be identified by various techniques ( [START_REF] Fisher | Statistical analysis of spherical data[END_REF], [START_REF] Marcotte | Automatic joint set clustering using a mixture of bivariate normal distributions[END_REF]). More generally, a cluster algorithm should be employed to determine fracture sets within which it is reasonable to assume independency between orientation and size distributions. Note that when the well radius is small compared to the fracture size, the size distribution itself is not needed to apply the lineic model. However, the latter is explicitly needed when the fracture size is of the same order of magnitude as the well radius. Unlike the orientation, the determination of the fracture size distribution remains a very difficult task. Indeed, well data are generally too scarce to infer precise information about the fracture 3D size and the only available methods to relate the latter to the trace distribution concern traces on plane outcrops (e.g. [START_REF] Berkowitz | Stereological analysis of fracture network structure in geological formations[END_REF], [START_REF] Priest | Determination of discontinuity size distributions from scanline data[END_REF], [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF] and [START_REF] Tonon | Closed-form and numerical solutions for the probability distribution function of fracture diameters[END_REF] for circular fractures, [START_REF] Zhang | Stereological relationship between trace length and size distribution of elliptical discontinuities[END_REF] for elliptic ones and [START_REF] Jimenez-Rodriguez | Inference of discontinuity trace length distributions using statistical graphical models[END_REF] for a method to estimate the trace distribution itself). Even if these methods can be extended to the case of a cylindrical well surface if the traces are small enough compared to the well radius, they can not be used if this last condition is not satisfied, which is often the case concerning pervasive fractures involved in fluid flows. Another attempt to characterize the 3D size distribution from well data developed in [START_REF] Özkaya | Fracture length estimation from borehole image logs[END_REF] is based on the assumption of a known (e.g. linear) relationship between the size of circular fractures and observable quantities such as the aperture or the spacing. Nevertheless, the issue of finding a robust method providing as precise information as possible about the 3D size distribution from well data without any assumed correlation with other variables may still need to be studied. Finally let us add that, among the basic assumptions allowing to derive the results of this paper, the Poisson spatial distribution of fracture centers deserves to be released. Indeed, although this distribution has been identified from a large number of field study ( [START_REF] Priest | Discontinuity spacings in rock[END_REF], [START_REF] Rives | Joint spacing: analogue and numerical simulations[END_REF], [START_REF] Sisavath | Geometry, percolation and transport properties of fracture networks derived from line data[END_REF]), it is suspected to be limited to low density networks [START_REF] Rives | Joint spacing: analogue and numerical simulations[END_REF]. It would thus be very interesting to consider correlations between fracture center positions such as in fractal models [START_REF] Darcel | Stereological analysis of fractal fracture networks[END_REF] or fracture networks such that the spacing on a well path follows other laws than the exponential one (11) characteristic of a Poisson process (e.g. lognormal in [START_REF] Rives | Joint spacing: analogue and numerical simulations[END_REF]).

A Calculation of surfaces

It is recalled that R c denotes the well radius and  and B =  cos Λ with Λ ∈ [0, π/2] respectively the major and minor radii of the ellipse obtained by projection of the fracture onto the plane perpendicular to the well direction. The section 3 highlights the correction factors involved when sampling fractures on a cylinder. The determination of the correction factors is based on the calculation of the following surfaces: the surface σ of the domain to which must belong the center of the projected elliptic fracture to cut the well circle (projection of the cylinder boundary) [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF], the surface σ f of the domain to which must belong the center of the projected fracture to contain entirely the well circle [START_REF] Lisle | Detection of zones of abnormal strains in structures using gaussian curvature analysis[END_REF] and the surface σ 2 of the domain to which must belong the center of the projected fracture to leave two traces on the well circle [START_REF] Narr | Estimating average fracture spacing in subsurface rock[END_REF]. The point of view in these definitions is that of an ellipse moving in the reference frame of the well circle. As suggested in [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], it reveals convenient to adopt the point of view of a circle moving in the reference frame of a motionless ellipse. It follows that, for instance, σ can also be seen as the area of the domain to which the circle center must belong in the reference frame of the elliptic surface so that the latter cuts it. Indeed, the two domains involved in the two definitions of σ are identical but centered respectively on the circle center and the ellipse center. The same reasoning can hold for σ f and σ 2 . Eventually, the determination of the boundary of the domains involved here suggests considering extreme cases with the curves described by the center of the circle when the latter is tangent to the ellipse.

A.1 Curves defined by a circle rolling on an ellipse

Two cases can be considered: either the circle "rolls" on the external side of the ellipse boundary or it "rolls" on the internal side of the ellipse boundary. Let us introduce the following parameterization of the ellipse M (α), its corresponding outward normal vector N (α) and the curvature radius ρ(α) with α ∈ [-π, π]:

M (α) = Â    sin α cos Λ cos α    ; N (α) =     cos Λ sin α √ 1-sin 2 Λ sin 2 α cos α √ 1-sin 2 Λ sin 2 α     ; ρ(α) = 1 -sin 2 Λ sin 2 α 3/2 cos Λ Â (59)
The curvature radius reaches its minimum at ρ(π/2) = Â cos 2 Λ and its maximum at ρ(0) = Â/ cos Λ. The curves described by the center of a circle of radius R c rolling outside (m + ) or inside (m -) the ellipse are parameterized by:

m ǫ (α) =     Â + ǫ Rc cos Λ √ 1-sin 2 Λ sin 2 α sin α Â cos Λ + ǫ Rc √ 1-sin 2 Λ sin 2 α cos α     ; ǫ ≡ + or - (60)
where ǫ ≡ + for the outside curve and ǫ ≡ -for the inside one. It is worth remarking that the parametric curves defined by (60) admit the x and y axes as symmetry axes, which means that only α ∈ [0, π/2] should be studied. To calculate surfaces bounded by those curves or parts of them, it will prove useful to introduce

j ǫ (α) = 4 × - 1 2 det m ǫ (α), m ǫ′ (α) (61)
which rewrites, using (60)

j ǫ (α) = 2 cos Λ Â2 + 2 cos Λ 1 -sin 2 Λ sin 2 α R 2 c + ǫ 4 1 -sin 2 Λ sin 2 α - 2 sin 2 Λ cos 2α + sin 2 Λ sin 4 α 1 -sin 2 Λ sin 2 α 3/2 ÂRc (62)
Surfaces bounded by the curves (60) will be obtained thanks to the following primitive of (62):

J ǫ (α) = 2α cos Λ Â2 + 2 arctan (cos Λ tan α) R 2 c + ǫ 4E (α, sin Λ) - sin 2 Λ sin 2α 1 -sin 2 Λ sin 2 α ÂRc ( 63
)
where

E (α, k) = α 0 1 -k 2 sin 2 θ dθ and E (k) = E (π/2, k) (64)
are respectively the incomplete and complete elliptic integrals of the second kind. The need to take care about the consistency between the sign convention adopted in (61) and the orientation of the domain boundary in surface calculations is worth being highlighted. On the one hand, the curve described by the center of any circle rolling outside the ellipse m + ([-π, π]) (60) is closed, regular, does not exhibit loops and is oriented in the clockwise direction when α increases. Therefore, considering the negative sign in the definition (61), the surface of the domain surrounded by this curve is simply provided by J + (π/2) in (63). On the other hand, the trajectory of the center of a circle rolling inside the ellipse may be described by different types of curves depending on the circle radius (see Fig. 7).

A.2 Calculation of σ

From now on, the reference frame is that centered on the circle center. In particular, the curves defined in the last paragraph A.1 are centered on the circle center and correspond to the motion of the center of the ellipse tangent to the circle. Recalling a result of [START_REF] Mauldon | Fracture sampling on a cylinder: From scanlines to boreholes and tunnels[END_REF], three cases have to be distinguished to calculate the surface of the domain to which the center of the ellipse must belong to intersect the circle. Finally the expressions of σ f are summarized in [START_REF] Lisle | Detection of zones of abnormal strains in structures using gaussian curvature analysis[END_REF]. Finally the expressions of σ 2 are summarized in [START_REF] Narr | Estimating average fracture spacing in subsurface rock[END_REF].

Figure 1 :

 1 Figure 1: Volume around the well path

10 E 90 Figure 3 :

 10903 Figure 3: Monte Carlo simulations compared to expectations and 1 st and 9 th deciles of P 32 given the trace number on the 1D well

Figure 4 :

 4 Figure 4: Euler angles of a fracture ellipse

Figure 5 :

 5 Figure 5: P.d.f. p P32|NL taking into account the double trace correction (solid line) or considering that all the traces come from different fractures (dashed line) and normalized histogram of P 32 for which discrete fracture networks have been built satisfying the same number of traces (here 20) along L = 20 m

M 1 model M 2 model M 3 Figure 6 :

 1236 Figure6: Monte Carlo simulations compared to expectations and 1 st and 9 th deciles of P 32 given the trace number on the well -the solid line (M 1) takes into account the double traces -the dashed lines (M 2) ignores the double traces -the dotted line (M 3) corresponds to the lineic result[START_REF] Dershowitz | Rock Joint System[END_REF]

  (a) If R c ≤  cos 2 Λ (min. curvature radius), the dashed curve of Fig. 7(a) described by the center of the circle is closed, regular, does not exhibit loops and is oriented in the clockwise direction. It therefore bounds a domain of area J -(π/2). (b) If R c ≥  cos 2 Λ (min. curvature and R c ≤  cos Λ (min. radius), the dashed curve of Fig. 7(b)shows singular points and loops. Invoking again the symmetries allowing to focus on α ∈ [0, π/2], it comes out that the point where the curve intersects itself, obtained for α = ξ, belongs to the x axis or equivalently m - y (ξ) = 0 in (60). This condition implies:ξ = arcsin Â2 cos 2 Λ -R 2 c  cos Λ sin Λ(65)At this point α = ξ, the circle is tangent to the ellipse in two symmetrical points with respect to the x axis M (ξ) and M (π -ξ).(c) If R c ≥  cos Λ (min. radius) and R c ≤  (max. radius), the dashed curve of Fig.7(c) shows singular points but no loops. Furthermore it is oriented in the counter-clockwise direction and thus defines a domain of area -J -(π/2).(d) If R c ≥  (max. radius) and R c ≤ Â/ cos Λ (max. curvature radius), the dashed curve of Fig.7(d)shows singular points and loops. The situation is similar to (b) except that the intersection point belongs to the y axis. The angle η ∈ [0, π/2] for which this point is reached is then defined by m - x (η) = 0, which implies:η = arcsin Â2 -R 2 c cos 2 Λ Â sin Λ(66)At this point α = η, the circle is tangent to the ellipse in two symmetrical points with respect to the y axis M (η) and M (-η).(e) If R c ≥ Â/ cos Λ (max. curvature radius), the dashed curve of Fig.7(e) is closed, regular without loops and oriented in the clockwise direction. It therefore defines a domain of area J -(π/2).

  (a) If  ≥ R c , the ellipse cuts the circle if and only if the ellipse center belongs to the domain bounded by the outside curve m + ([-π, π]) (60) (see Fig. 8(a)). Consequently, σ = J + (π/2). (b) If  ≤ R c cos Λ, the ellipse center must belong to the domain bounded by m + ([-π, π]) but must be outside the domain bounded by m -([-π, π]) (see Fig. 8(b)). Indeed, when the ellipse center belongs to the latter, the ellipse is entirely contained in the circle but does not cut it. The surface is then σ = J + (π/2) -J -(π/2).

Figure 8 :

 8 Figure 8: Domains of ellipse centers leading to a non empty intersection

(a) Â cos 2 ΛFigure 9 :A. 4 Calculation of σ 2 If

 2942 Figure 9: Domains of ellipse centers leading to full traces ( Â cos Λ > R c )

Figure 10 :

 10 Figure 10: Domains of ellipse centers leading to double traces ( Â cos 2 Λ < R c < Â/ cos Λ)

with η given by (66). Indeed, the two other domains correspond to ellipses which cut the circle along two traces and must be taken into account in σ. Hence, the surface is

Finally those three results are summarized in [START_REF] La Pointe | Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations[END_REF].

A.3 Calculation of σ f

A full intersection means that the circle is entirely contained inside the elliptic domain. The first necessary condition is obviously that the minimum radius of the ellipse  cos Λ be greater than R c otherwise σ f = 0. Considering now  cos Λ > R c , the external boundary of the domain to which the ellipse center must belong is then the curve described by m -(α) (60). Indeed, the extreme case is obtained by situations where the circle and the ellipse are tangent but it is also necessary to keep on satisfying the fact that circle should stay inside the ellipse. Two cases can arise: