Non-parametric three-dimensional discrete fracture network modelling from fracture mapping data

Rock masses usually include discontinuities in the form of fractures or joints, which in general are critical in determining the mechanical or flow properties of the rock mass at the engineering scale. Recent advances in computing power have made it possible to represent these fractures explicitly in a rock mass model. However, in practice, three-dimensional fracture systems are impossible to measure, and the best available alternative is either to map fractures in drill cores (one-dimensional) or on exposed rock surfaces (two-dimensional). The construction of a three-dimensional fracture system that is representative of reality then becomes a very challenging inverse problem, as this is essentially a relationship of one reality corresponding to many possibilities. Most existing approaches start with assumptions about the distribution of measured features, perform a series of bias corrections and then use some further assumptions and stereological analysis to establish the three-dimensional fracture model. The usual practice is to focus on particular aspects (parameters) of the three-dimensional fracture model which can easily produce a biased representation of the real fracture system. In this paper, a non-parametric systematic optimisation framework is described which can help construct a more realistic three-dimensional fracture model that best matches the statistics of the fracture mapping data. No bias correction is necessary in this approach provided the sampling of the fracture model is consistent with that used in the data collection. The framework can also be used to assess the uncertainty of the fracture model, which is an important, but often neglected, component of modelling in existing approaches. In addition, this non-parametric approach is also flexible and can deal with data from different sources (e.g., different sampling planes at different orientations, mixtures of plane and borehole sampling). This is a significant advantage over existing approaches as different data sources often display different statistical characteristics even for the same fracture set. The Stripa Mine dataset is used to demonstrate the application of the proposed method.


Introduction
In general, a rock mass consists of intact rock and fractures/joints.The intact rock component refers to the cemented mineral grains plus the voids between grains and is normally modelled as a continuous medium.In the context of this work, potential micro-cracks existing within this "intact" rock structure are not explicitly identified and are only treated as part of the "intact" structure.Fractures on the other hand represent discontinuities between intact rocks where separation boundaries are clearly observable with the naked eye.In the context of the mechanical behaviour of rock masses, these discontinuities are usually the weakest links in the system and can cause serious stability issues, while in the context of the flow behaviour of rock masses, these discontinuities are usually the major conduits that conduct flow through the system as, in general, the permeability of fractures could be orders of magnitude higher than that of the rock matrix (intact rock), particularly in hard rock settings.In other words, fractures may occupy only a small fraction of the total volume in a rock mass, but they play a decisive role in determining the mechanical and flow behaviours of rock masses.
It is therefore important that the fracture system is properly represented in any rock mass model used to assess its engineering behaviour.In practice, a three-dimensional (3D) fracture system cannot be measured directly and only a partial image of the system is observable, either in the form of drill-core surface images or two-dimensional exposed rock faces.Although various scanning techniques and computed tomography imaging methods are commonly used to analyse rock samples, the scale is, in general, too small to have much significance for engineering applications.3D fracture mapping based on cut slices through a large block of rock have been conducted (Dowd et al. 2015) but this is a costly process that is extremely time-consuming and is difficult to adapt to wide-spread applications.Other types of data sets of fracture systems may also be available from specific applications (e.g., seismic event points from engineered geothermal reservoirs, see, for example, Xu et al., 2013, Xu andDowd, 2014).
Data from fracture mapping of drill cores are essentially one-dimensional (1D).Orientations of fractures (dip directions and dip angles) can be measured using oriented cores.Fracture spacing can also be readily measured along drill cores.This data set, however, is seriously biased due to the orientation of the drill holes although there are techniques for correcting the bias in the statistics of fracture spacing (e.g., Terzaghi, 1964).Data from fracture mapping of exposed rock faces are two-dimensional (2D) and can provide much more information than 1D data.Scanline surveys can be applied to measure fracture spacing (e.g., Priest and Hudson, 1981) and both scanline and window surveys can be used to measure fracture trace length (Song andLee, 2001, Priest, 2004), which is the most important signature of 3D fracture sizes.There is also bias in 2D datasets due to the orientation of the sampling faces in relation to the inner fracture system.Fracture trace length measurement may also be biased due to truncation (ignoring small fractures) or censoring (one end or both ends of a fracture trace are not visible) because of the limited size of the sampling face (hence the survey window).Again, tools have been developed to make corrections in the derived statistics based on certain assumptions (e.g., Laslett, 1982).Both rectangular and circular sampling windows have been used in published research (Zhang and Einstein, 2000) but the latter option naturally avoids the bias issue due to the orientation of the sampling window.In general, fracture spacing is used to estimate fracture density and trace length is used to estimate fracture size, although fracture size estimation from drill-hole image logs has also been attempted (Özkaya, 2003).
Fracture mapping data provide the basic inputs for the derivation of the parameters of the 3D fracture system.The first step of the most common approach is to classify mapped fractures into different clusters, representing different fracture sets, based on the similarity in fracture orientations (Priest, 1993).The classification is based on the belief that fractures formed by the same geological or mechanical process will, in general, have similar orientation characteristics.Lower hemispherical fracture pole projection is the common tool used in this analysis.Although various automatic or semi-automatic clustering algorithms have been applied with various degrees of success (Mahtab and Yegulalp, 1982, Hammah and Curran 1998, Zhan et al., 2017), manual clustering based on visual inspections remains an effective way for fracture set classification.After identification of fracture sets, 3D fracture size and fracture density parameters are derived on a set by set basis, although in some cases correlations between sets can also be considered (Chiles, 1988, Dowd et al. 2007, Long and Billaux, 1987).Current approaches usually require assumptions about the statistical distributions of the measured features (e.g., fracture trace length and fracture spacing), 3D fracture size and fracture locations so that their parameters can be estimated.Bias corrections are commonly performed first on the distribution parameters of the measured features before they are used to derive the 3D fracture parameters.Additional assumptions to simplify the 3D fracture system are also needed (e.g., uniform distribution of fracture locations) in order to establish the relationship between the 3D system and the 2D or 1D mapping data based on stereological analysis.These assumptions are necessary to reach a solution, but they do not necessarily reflect accurately the reality of the mapped fracture system.
This paper describes a non-parametric optimisation process that can be used to estimate the parameters of the 3D fracture system based on the statistics of 2D or 1D fracture mappings.The process is essentially a minimum c 2 method, commonly used in statistics to maximise the goodness-of-fit of target distributions.
No assumption about the distributions of the measured features from the 1D or 2D fracture mappings is required and there is no correction of the sampling biases.In other words, the biased histogram, based on the limited size of the scanline or sampling window, is used directly in the parameter fitting process.The method can also accommodate non-homogeneous spatial distributions of fractures, i.e., the common assumption in existing approaches of a uniform random distribution of fractures in 3D space is no longer necessary.In practical applications, different types of fracture mappings are commonly employed resulting in a dataset comprising different sources (e.g., different mapping done on sampling planes with different orientations, a combination of mappings done on rock faces and rock cores), see for example, Rouleau and Gale (1985a).A common problem encountered when using different sources is that the relative date may produce different statistics even for the same fracture set, which, for parametric approaches, will pose a serious reconciliation problem.Kulatilake et al. (1993) noted that there are significant differences between the estimated statistical parameters based on the fracture mapping data of the floor and the walls of the Stripa ventilation drift, which are likely to be caused by orientation bias.The non-parametric approach described in this paper is flexible in the sense that it can incorporate different biased sampling statistics into the optimisation process so that a global optimal solution can be found.The fracture mapping dataset of the Stripa Mine ventilation drift are used in this study to demonstrate the application of the proposed non-parametric method.
In the following sections, existing approaches to the derivation of parameters of a 3D fracture model are briefly reviewed, followed by a description of the minimum c 2 method.The data set used in this study is then introduced, together with a summary of target statistics.The proposed method is then applied to the case study data set and the parameters of the 3D fracture model are derived.The final model is generated and the overall statistics are then analysed and compared with those obtained from the mapping data.

Fracture size
The assumptions are: 1).fractures can be represented as circular discs with diameter D; 2). the centres of fractures are uniformly randomly distributed in 3D space and the volume density follows a Poisson distribution; 3).fractures in the same set are parallel; 4).there is no correlation between fracture size and fracture location; and 5). the sampling plane is infinite and thus there is no censoring bias.
Based on these assumptions, an analytical relationship between the statistical distribution of D, g(D), and the statistical distribution of the length of the interaction traces between fractures and the sampling plane, f(l), can be established based on stereological analysis (Warburton, 1980, Zhang andEinstein, 2000): (1) where is the mean of D. In practice, to make use of Eq. ( 1), it is important to correct f(l) to account for sampling biases.These biases mainly include (Zhang and Einstein, 2000;Einstein and Baecher, 1983): 1).censoring bias due to the limited size of the sampling window -long fracture traces may have one end invisible (half-censored) or both ends invisible (fully-censored); 2).truncation bias -small traces are not mapped; 3).orientation bias -relative orientations between sampling plane and fractures cause different fractures that have different probabilities of being intersected by the sampling plane; and 4).size bias -larger fractures have a higher probability of being sampled.The first two biases can easily be accounted for using statistical approaches and there is no shortage of proposed methods depending on the assumptions involved (see, for examples, Pahl, 1981, Priest and Hudson, 1981, Laslett, 1982, Kulatilake and Wu, 1984, Villaescusa and Brown, 1992, Zhang and Einstein, 2000, Song and Lee, 2001).Orientation bias is difficult to eliminate in practice as normally only a few sampling planes are available.However, at the very least, a reconciliation of the estimated parameters is necessary if the estimations are to be done on different sampling planes.Size bias is normally not incorporated directly in the estimation.Note that to make the bias correction for f(l), the type of distribution has to be assumed for f(l) and the bias correction will then focus on the correction of the parameters of the distribution, such as its mean and standard deviation.Types of distributions for f(l) commonly used in practice are exponential, lognormal and gamma distributions (Villaescusa and Brown, 1992, Zhang and Ein-stein, 2000, Tonon and Chen, 2007).Note that the trace length distribution f(l) can also be estimated from the scanline sampling.In this case, the relationship between f(l) and g(D) that is used is slightly different to that in Eq. (1) (Warburton, 1980).
Once the "true" f(l) is obtained, Eq. ( 1) can be used to derive the distributions for g(D).For analytical solutions, the type of distribution for g(D) must first be assumed.The three commonly used types of distribution for trace length are also commonly used for g (D).Based on the assumed distribution, the relationships between the parameters of g(D) and f(l), such as mean and standard deviation, can then be established.Examples can be found in Zhang and Einstein (2000), Tonon and Chen (2007) and Hekmatnejad et al. (2018) who used a distribution-free approach to solve the problem.If g(D) is not assumed, numerical solution techniques can be used to estimate the distribution directly (Kulatilake and Wu, 1986), which requires an iteration procedure to find the best solution.

Fracture intensity
Fracture intensity can be expressed using the persistence measure P ij introduced by Dershowitz (Dershowitz, 1984, Dershowitz andHerda, 1992) where i represents the dimension of the measurement region (1, 2 or 3) and j represents the dimension of the measured feature of the fracture (0, 1, 2 or 3).For example, P32 is the average fracture area per unit volume which is an important measure related to the degree of fracturing (fracture intensity) within the rock mass.A Discrete Fracture Network (DFN) model can be generated by matching different fracture intensities but the simplest implementation is based on P20 for 2D and P30 for 3D applications.P30 (P20) is the average number of fractures per unit volume (area), which is also commonly referred to as fracture density.This definition of fracture density is based on the convenience of representing the location of each fracture by a point on the fracture (normally the centre point of the fracture).The fracture distribution is then modelled as a point process and the DFN modelling is a marked point process as other fracture properties, such as fracture size and orientation, are treated as random marks associated with the points.The most commonly used point process model is the Poisson model where points (fracture locations) are uniformly (randomly) distributed in space, even though fracture distributions are generally non-homogeneous (Xu and Dowd, 2010).Note that P10 of a fracture set is the inverse of the fracture spacing in the normal (or average normal) direction of the fracture plane.
Fracture spacing is a common statistical measure obtained from scanline sampling along the drill cores or on exposed rock faces.If the scanline does not coincide with the normal direction of the fractures, a simple factor of cos a can be applied to correct the orientation bias, where a is the acute angle between the scanline and the normal direction of the fractures (Terzaghi, 1964, Fouché andDiebolt, 2004).Common types of distributions for fracture spacing are exponential, lognormal and gamma distributions, where exponential distribution arises naturally if fracture are uniformly distributed in space (Priest and Hudson, 1981).
In this work, fracture volumetric density, P30, is required for the DFN modelling.Oda (1982) established the relationship between P30 and P10 along a scanline as: (2) where r = D/2 is the radius of circular fracture disc, n is the unit vector representing the normal of the fracture plane, i is the unit vector representing the direction of the scanline and E[.] is the expected value.This relationship is used by researchers (e.g., Kulatilake et al, 1993) to estimate P30, i.e.: (3) P30= 4P10 π .E [D 2 ] .E[|n .i|] Kulatilake et al. (1993) also established a simple relationship to estimate P30 from P20, as given by: where a is the angle between the sampling plane and the fracture plane.In these approaches, a Poisson distribution of fracture location is assumed.Sampling bias due to the limited size of the scanline or sampling window must be considered when using these relationships.Note that the estimations discussed above are only for the expected density values.Other distribution characteristics are not considered here.

A general non-parametric estimation method
The c 2 test is commonly used to test the goodness-offit of a statistical distribution model.The c 2 statistic is calculated as: (5 where K is the number of mutually exclusive classes used to cover the range of values for the random variable, and are the observed and expected numbers (based on an assumed distribution) of values of the random variable that fall within class i. Statistic of Eq. ( 5) follows a c 2 distribution in the limit when the total number of observations becomes large.The statistic quantifies the total percentage difference between the observed and the tested distributions and therefore a smaller value means a better goodness-of-fit of the observation to the assumed distribution.
In this study, a version of the c 2 statistic is calculated for trace length and fracture spacing.These statistics are used to set up the primary objectives for the optimisation process to derive the parameters of the 3D fracture system.In this case, the tested distribution models (corresponding to the expected number of observations ) are in fact the target distributions obtained directly from the fracture mapping data.The observed distribution models (corresponding to the number of observations x i ) are the empirical distribution obtained by sampling the trial 3D fracture system (model).The modified c 2 statistics used in this application are calculated as: (6) where p i and q i are the relative frequencies of the fracture model and target distributions in class i.The calculated statistic is then used to quantify the goodness-of-fit of model to the target distribution.
The key differences between this non-parametric approach and the traditional approach described in the previous sections are that it is not necessary to assign a distribution type to the measured features (e.g., trace length or fracture spacing) and it is not necessary to correct for bias.The target distributions that are to be matched by the model stay in non-parametric forms.Of course, simulation is needed to sample the fracture model and this can be done in the same way as that used in the fracture mapping process.
The main characteristics of the 3D fracture system that need to be derived from the optimisation are the fracture size distribution g(D) (or g(r)) and the fracture density P30(X), which in general depends on location X (i.e., a non-homogeneous spatial distribution).Although in this study both the size distribution and the density are found for each fracture independently (i.e., assuming no correlation between fracture sets), the method can easily be extended to cover fracture set correlations such as the hierarchical model described in Lee et al. (1990) or the cluster model described in Xu and Dowd (2010).The flowchart of the proposed optimisation process is shown in Fig. 1.
For each set of g(D) and P30(X) models, M number of independent realisations are used to derive reliable statistics, including the range of variations which can later be used to quantify the uncertainty of the estimates.This is repeated for all sets of models and the set with the minimum c 2 statistics will be regarded as the optimal solution.A brute-force search scheme can always be used if the combined search space for the model parameters is limited.If the number of parameters is significantly high, smarter search techniques such as Markov chain simulation can be used (Xu et al. 2013).If the spatial distribution of fractures is homogeneous (i.e., constant P30), or if there is no correlation between fracture size and fracture location, the distribution of fracture trace length will be independent of fracture density (see also Eq. ( 1)).In this case, the optimisation process can focus first on g(D) for a set fracture density value.The optimal solution found for g(D) is then used to optimise the P30(X) model.However, fracture size and fracture density could be correlated.For example, there may be a systemat- ic change in terms of fracture density and fracture size against depth (Xu and Dowd, 2010).In this case, the parameter space defined by the range of [g(D), P30(X)] must be fully searched jointly to find the optimal solution.Note that although the fracture size is written as g(D), the proposed optimisation framework is completely general for different types of fracture geometries.Observations of both equidimensional fractures (i.e., approximately circular) and non-equidimensional fractures are reported in the literature (Zhang and Einstein, 2000).If a geometry other than a circle is assumed for the fractures, the same optimisation process shown in Fig. 1 can still be followed with some added parameters in the DFN model.
One major advantage of this non-parametric optimisation approach is that the sampling bias corrections for the target statistics are no longer necessary because both the target and model statistics are biased in the same way (truncation, censoring or orientation).This does not necessarily mean that biases no longer need to be considered in this model optimisation process.In fact, if the sampling window is limited in size (such as in the case study described later), censoring can become serious and therefore the censoring statistics of fracture traces should also be built into the decision process.These statistics include the proportions of fractures traces with both ends visible (T0), one end visible (T1) and both ends invisible (T2).In this case, to ensure that target and model statistics are comparable, the same sampling window must be used for the fracture mapping data and for sampling the DFN models.(Rouleau et al. 1981).Figura 2. Área experimental en Stripa Mine (Rouleau et al., 1981).

A case study -the Stripa Mine fracture dataset
The Stripa Mine is an old iron ore mine located approximately 145 km west of Stockholm, Sweden.It was the test site for a comprehensive research and development programme on nuclear waste storage in fractured granitic rocks in the 1970s and 80s (Rouleau et al. 1981, Rouleau and Gale, 1985a, 1985b).The experimental area at Stripa included several drifts in the crystalline rock 340 m below the surface, as shown in Fig. 2. The ventilation drift, approximately 45 m long, located at the north of the test area, was excavated mainly for the fracture hydrology experimental programme.For this purpose, detailed fracture mapping was performed in the major section of the drift, including mappings of the floor and the two walls.A typical map of fracture traces is shown in Fig. 3.All fractures are numbered, and their strike and dip directions are recorded together with their length rounded to the nearest metre, aperture, infill materials and fracture roughness.A truncation threshold of 0.5 m was used and only fracture traces greater than the threshold were mapped.The measurements are coded in a text file available in Rouleau et al. (1981).Part of this dataset (a section of 10m) was used in Kulatilake et al. (1993) to demonstrate their technique to derive the fracture model for the mine site using the traditional parametric approach discussed above.
A total of 279 fracture traces were mapped on the East wall, which is 36 m long; 338 fracture traces were mapped on the West wall which is 56 m long; 371 fracture traces were mapped on the floor, which is 42 m long; and 38 fracture traces were mapped on the northern end face of the drift.Note that there are densely fractured zones (e.g., areas with dotted lines in Fig. 3) with fracture trace lengths < 0.5 m and fracture spacing ranging from less than 5cm to 10 cm in the area, which are not considered herein.Poles of the mapped fractures were plotted using lower hemispherical equal-area projection for fracture set classification, as shown in Fig. 4.
Clearly there are three major clusters, sets 1, 2 and 3. Set 4 is not significant, but it does show a certain degree of concentration.Set 5 is the background set which account for all fractures not assigned to the four clusters.Note that this fracture set classification is slightly different to that documented in Rouleau and Gale (1985) for the ventilation drift where sets 2 and 3 shown in Fig. 4 were combined into one set.Fractures not assigned to any of the clusters were ignored in their classification (i.e., no background set), although these fractures make up about half of all the fractures mapped (see Table 1).The features of the clusters shown in Fig. 4 were also further confirmed by incorporating an additional 1622 fractures mapped from boreholes drilled in the ventilation drift (Rouleau and Gale, 1985a).For the classification shown in Fig. 4, the orientation statistics are listed in Table 1.Once fracture sets are classified, fracture trace maps such as the one shown in Fig. 3 need to be digitized for different fracture sets.The trace length recorded in the original data file was rounded to the nearest metre and the location of the trace was recorded as the ID of areas where the fracture trace was first encountered.This information is not adequate to derive the statistics of trace length and fracture spacing.The digitised trace maps for each set of fractures are shown in Fig. 5.The length of the mapped east wall is set at 37 m, the floor at 43 m and the west wall at 57 m.The vertical height of the mapped wall or the width of the mapped floor displayed is 4 m.Note that the arched profile of the tunnel walls is ignored in this work.The estimated P20 for each case is also given.

Set
The fractures in fracture set 1 are sub-vertical with a mean dip angle of 71 o , dipping approximately at an azimuth of 18 o (left boundary of Fig. 5 is the north face of the ventilation drift).This set of fractures shows up reasonably well on all three mapping planes with the fracture density (P20) of the floor slightly higher than that of the walls.The fractures can be considered approximately evenly distributed.The scanlines to be used will be oriented along the drift direction, i.e., parallel to the axis of the tunnel, which will ensure that most fractures will be intersected.Note that in the method presented here, it is not absolutely necessary for the scanlines to be perpendicular to the mean fracture plane of the fracture set.No correction of the scanline orientation bias is necessary as long as the scanlines can intersect fractures at a reasonable angle.The sizes of the traces shown on the three mapping planes are very similar which suggests that an equidimensional assumption (i.e., circular disc) is acceptable for this set of fractures.
Fracture set 2 comprises sub-horizontal fractures with a mean dip angle of 20 o , dipping approximately towards the East (going into the page on the East wall images).The fracture traces on the two walls are consistent in terms of the orientation of the fracture Figure 6.Trace length distributions (horizontal axis is trace length l in metres and vertical axis is relative frequency).Figura 6. Distribuciones de longitud de trazas (el eje horizontal es la longitud de la traza l en metros y el eje vertical es la frecuencia relativa).
set.For the floor, the mapping plane orientation bias appears to be significant in this case as the characteristics of the traces are quite different to those shown on the two walls.In terms of the spatial distribution of fractures, set 2 has a higher degree of non-homogeneity than set 1.However, a uniformly distributed model is still used in this exercise for reasons of simplicity.Note that towards the south of the West wall (right side of West wall images), there is a cluster of fracture traces.This cluster may cause an overestimation of the fracture density for this fracture set (see later discussion).The orientation of the scanlines will again be parallel to the axis of the tunnel to ensure the maximum number of intersections, even though approximately vertical scanlines should be used for this fracture set to minimise scanline orientation sampling bias.As the optimisation technique described in this paper does not depend on the statistics corrected for sampling biases, it is more important to use scanlines that have more fracture intersections so that more reliable statistics can be obtained, even if they are biased.An equidimensional assumption will again be used for fracture size.
Fracture set 3 comprises sub-horizontal fractures with a mean dip angle of 19.5 o , dipping approximately towards the south (right side of the images shown in Fig. 5).Although the fracture traces shown on the three mapping planes are consistent with the orientation of this set of fractures, a significant issue is that the fracture density on the east wall is approximately four times higher than the fracture densities of the other two mapping planes.This is inconsistent with the assumption of a random distribution and suggests that a non-homogeneous fracture distribution model is more appropriate for this fracture set.However, as there is not sufficient fracture mapping to define the variation Figura 7. Fracture spacing distributions (horizontal axis is fracture spacing in metres and vertical axis is relative frequency).

Figura 7. Distribuciones de espaciado de fracturas (el eje horizontal es el espaciado de fracturas en metros y el eje vertical es la frecuencia relativa).
of fracture distributions in the East-West direction, an average fracture density is used in this exercise.
Fracture set 4 has a dip angle of 60 o , dipping approximately towards the West (coming out of the page on the East wall images).Fracture traces shown on the three mapping planes are consistent with the orientation of the fracture set.The number of fractures for this set of fractures is small, particularly for the East wall.
Fracture set 5 is the background set with assumed purely random orientations.The fracture density values of the East wall and the floor are similar but that for the West wall is much lower.This, together with the features of fracture set 3, suggest that the East wall is much more fractured than the West wall (overall P20 = 1.8851, 2.1570, 1.4825 for East wall, floor and West wall) and therefore a non-homogeneous fracture distribution model is more appropriate.However, as mentioned above, two fracture mapping planes in the East-West direction are not sufficient to define a non-homogeneous model and therefore an overall average homogeneous model is used in this exercise.

Target Statistics
The target statistics for the model to match in our approach are raw sampling statistics without any bias correction.As different sampling planes with different orientations will, in general, produce different sampling statistics, two sets of target statistics will be derived in this case: one from the fracture mappings of the two walls (vertical sampling plane) and the other from the floor (horizontal sampling plane).Statistics of two fracture characteristics are required, one for the    trace length used to optimise the fracture size parameters and the other for the fracture spacing to be used to optimise the 3D fracture density.For trace length, the rectangular sampling window is defined as the size of the fracture mapping window, i.e., 4 m x 37 m for the East Wall, 4 m x 43 m for the floor and 4 m x 57 m for the West Wall.For fracture spacing, four scanlines parallel to the axis of the tunnel are used for fracture sets 1 and 5 for each mapping face, located at 0.5 m, 1.5 m, 2.5 m and 3.5 m along the width of the sampling window.For fracture sets 2, 3 and 4, where four scanlines are not sufficient to derive the required statistics (too few intersections with fracture traces), the number of scanlines is increased systematically until reliable statistics on fracture spacing are obtained.The target statistics for trace length are shown in Fig. 6 and the target statistics for fracture spacing are shown in Fig. 7.
Clearly, the trace length distributions are positively skewed, which is consistent with many observations in practice (Zhang andEinstein, 2000, Hudson andPriest, 1983).Most published datasets suggest the type of distribution as exponential, lognormal or gam-ma, although fitting a model is not entirely necessary in the non-parametric approach described in this paper.For the Stripa dataset, the trace lengths from the two walls tend to follow a lognormal distribution for all fracture sets, while the trace length obtained from the floor mapping suggests an exponential distribution.As shown in Fig. 6, a much greater proportion of small fracture traces are observed on the floor (horizontal plane) than on the two walls (vertical plane).This suggests that the 3D fractures in this mine site may not follow exactly the equidimensional assumption, at least for some of the fracture sets.However, as the overall trends are still largely similar, for simplicity the equidimensional assumption is used in this work to demonstrate the application of the proposed method.In terms of the target statistics, the average of the two walls will be used for the vertical sampling plane and that of the floor will be used for the horizontal sampling plane.Fig. 8a shows a comparison of the average target statistics of the trace lengths of different fractures sets.From the figure, it is reasonable to assume that the statistical characteristics of the frac-Figura 9. c 2 values of the fracture sets at different l values (horizontal axis is l and vertical axis is c 2 ).Figura 9. Valores c 2 de los conjuntos de fracturas para diferentes valores l (el eje horizontal es l y el eje vertical es c 2 ).ture lengths of different sets of fractures are largely the same.This suggests that the fracture sizes of different fracture sets may, in this case, follow the same or a similar distribution with similar parameters.This observation will be further explored in the parameter optimisation process discussed later.
The fracture spacing distributions are positively skewed, which is as expected from extensive published observations (Priest and Hudson, 1981).The types of distributions commonly used in published research include exponential, lognormal and gamma, with the exponential as the most popular as it is the de facto type of distribution if the fracture distribution in space is purely random.For the Stripa data set, fracture spacing for Sets 1, 2 and 5 are clearly exponential, but not Sets 3 and 4.This is consistent with the discussion on the fracture distributions in Fig. 5, although there is a smaller number of fractures in Set 4. For the target statistics to be used in the parameter optimisation process, the average of the two walls will be used for the vertical sampling plane and the floor will be used for the horizontal sampling plane.A comparison of the average target statistics of fracture spacing for the different fracture sets is given in Fig. 8b, which clearly shows significant differences.Distributions for fracture sets 2, 3 and 4 are similar, although fracture set 3 shows a significant upper tail for spacing greater than 20 m.Sets 2 and 4 are very similar although set 2 has a much higher proportion of smaller fracture spacing.Fracture set 5 has the smallest fracture spacing, followed by fracture set 1. From Fig. 8b, it is reasonable to expect that the fracture density for fracture set 5 will be the highest, followed by set 1 and then set 2, 3 and 4, under the assumption of a random distribution of fracture orientations.As this assumption is certainly not true, these expectations may not correctly reflect the actual fracture density values.Sets 2, 3 and 4 should have similar density values although the density of set 3 should be slightly lower -without considering orientation bias.
The additional statistics that will be used in the parameter optimisation process are the fracture censoring statistics listed in Table 2, which show that there are almost no long fracture traces for which both ends are censored.The majority of traces are fully contained within the defined sampling window (overall average

> 80%
).This suggests that fractures in the Stripa Mine are, in general, small in size with the average much less than 4 m.Note that the censoring statistics are only used when there are difficulties in making decisions when c 2 statistics are used alone.

Fracture Size
As discussed above, an equidimensional fracture model will be used to demonstrate the parameter optimisation process to derive the parameters for the distribution function g(r), where the radius of the fracture disc instead of the diameter is modelled.As a simple example to demonstrate the technique, an exponential model is used for which only one parameter (l) is required, i.e., , where the mean fracture size is 1/l.For fracture density, under the assumption of a random distribution, the relationship between f(l) and g(r) does not depend on P30.Therefore, P30 is set to 1.0 in the optimisation for the fracture size distribution parameters.
Fig. 9 shows the variation of c 2 for the fracture sets at different l values, where c 2 is calculated for the trace length target statistics shown in Fig. 6 and by sampling 3D DFN models.The volume used for the 3D DFN simulation is 37 m x 6 m x 6 m, where the enlarged cross-sectional area (6 m x 6 m) is used to reduce the potential boundary effect.The sampling plane is 37 m x 4 m, either oriented vertically simulating the wall sampling or horizontally simulating the floor sampling.A total of 40 independent simulations were generated for each set of parameters to quantify the 95% confidence interval of the derived statistics.
The l values for fracture set 1 optimised for the wall and floor target statistics are both at l = 2.2.This set of fractures intersect the axis of the tunnel at a high angle and, as can be seen from Fig. 5, fracture traces for this set of fractures were well captured in both the wall and floor mappings and therefore the optimisation from both target statistics is expected to yield similar answers.For fracture set 2, as can be seen from Fig. 5, the wall and floor target statistics are significantly different as are the optimised results shown in Fig. 9.As this set of fractures is sub-horizontal, some differences are expected due to the sampling plane orientations, but the amount of the difference is unexpected.Based on the wall target statistics, l = 2.2, but based on the floor target statistics l = 8.8, which is clearly unrealistic in this case.For the censoring statistics, at l = 2.2, (T0, T1) = (93.8%,6.2%) and (89%, 11%) for wall and floor sampling respectively, compared with the target statistics of (93.4%, 6.6%) and (70.4%, 29.6%) respectively.For l = 8.8, (T0, T1) = (98%, 2%) and (99%, 1%) for wall and floor sampling, which indicates a much worse fitting in terms of censoring statistics.As discussed above, the similarity of the distributions of trace lengths suggests similar statistical characteristics of fracture sizes for different fracture sets, therefore l = 2.2 based on wall mapping was chosen for this fracture set.For fracture set 3, the same problem is encountered due to significantly different target statistics derived from wall and floor mappings and the final optimised value is the one derived from the target statistics of wall mapping, i.e., l = 2.2.Fracture set 4 has only a small number of fractures.The optimisation based on the wall mapping gives l = 2.2 and that based on the floor mapping gives l = 2.0.Again, the value of l = 2.2 is chosen.This parameter gives the censor statistics of (T0, T1) = (85%, 15%) and (91%, 9%) for wall and floor sampling planes, compared with the corresponding target censoring statistics of (82.5%, 17.5%) and (82.2%, 13.8%) respectively.For fracture set 5, the optimised l based on wall statistics is 3.4 and that based on floor statistics is 4.8.As this background set of fractures has a random distribution in terms of fracture orientation, the statistics from the wall and the floor mappings can be combined.The optimisation using the combined statistics gives l = 2.8, which will be the chosen l value for this fracture set.In terms of censoring statistics, this l value gives (T0, T1) = (78%, 22%) for vertical and (74%, 26%) for horizontal plane sampling, which compares favourably with the corresponding target statistics of (86%, 14%) for wall and (72.6%, 27.4%) for floor mapping.rises the final optimised parameter l for all five fracture sets.These optimised parameters are then used in subsequent optimisation for fracture spacing.

Fracture density
Similar procedures (cf.Fig. 1) will be followed to derive the optimal values for the fracture density for each fracture set, based on the fracture spacing target statistics shown in Fig. 7.
For fracture set 1, the optimised P30 value based on the fracture spacing for the wall mapping is 0.5 and that based on the floor mapping is 0.9.As this set of fractures dips approximately 70 o with a dip direction of approximately 18 o from North, the vertical and horizontal sampling planes will intersect the fractures at approximately the same angle.Therefore, similar features in terms of fracture intersections should be expected, which are evident in the patterns shown in Fig. 5.For this reason, the target statistics from wall and floor mappings can be combined to derive the optimal value for P30 and the result is P30 = 0.6 as shown in Fig. 10.For fracture set 2, the optimised P30 from wall and floor mappings is between 1.0 to 1.8.This set of fractures is sub-horizontal and significantly different features are expected from vertical and horizontal plane sampling due to the orientation bias.The target statistics, therefore, cannot be combined for the optimisation.After close inspection the difference between the model and target statistics, it was decided to use P30 = 1.0 since after that value there is no significant change in the goodness of fit.For fracture set 3, P30 was optimised based on the target statistics of 0.4 for wall mapping and 1.1 for floor mapping.Again, on close inspection of the differences between the target and model statistics, the compromise choice was P30 = 0.6.For fracture set 4, the optimised P30 values from the target statistics of both the wall and floor mappings are almost identical with the best choice being P30 = 0.8.For fracture set 5, the P30 based on wall mapping is 2.2 and that based on floor mapping is 2.6.As the orientations of this set of fractures follow a random distribution, both statistics can be combined.The optimised P30 based on the combined statistics is P30 = 2.6.The final optimised P30 values are listed in Table 3.Some of the optimised P30 values seem to be high, e.g., that for set 5. This is mainly caused by the large l value of the corresponding fracture set.For fracture set 5, l = 2.8 means that their average fracture radius is 0.36 m (or diameter 0.72 m).As a truncation threshold of 0.5 m was applied in fracture mapping, only 50% of the fractures in this set (with r > 0.25m) will have any chance of being recorded in the fracture trace maps.For this reason, a large number of fractures will have to be generated in order to have sufficient fracture intersections (with trace length > 0.5 m) to match that observed in the trace maps.Hence, P30 will be high.For the other fracture sets, the percentages of generated fractures that have sizes less than 0.5 m, and thus have no chance of being recorded, are 14% for sets 1 Figura 13.Comparación entre el modelo y estadísticos de los datos de longitud de la traza (parte superior) y espaciado de fracturas (parte inferior) de todos los conjuntos de fracturas (el eje horizontal es la longitud de la traza (parte superior) y espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa.and 3, 22% for set 2 and 18% for set 4. As there was no report of any significant number of small-size fractures in the drift, this is perhaps an indication that the exponential fracture size distribution is not a suitable assumption for this dataset.The Stripa fracture mapping work report (Rouleau et al. 1981) does mention that there is a large number of small-scale fractures (with trace length < 0.5 m and a spacing of 5 to 10 cm) observed on the walls and the floor, but these fractures occur in clusters termed fractured zones in the report.They are fundamentally different to small-scale fractures due to an assumed size distribution and are expected to be uniformly distributed in space based on the model used.For these reasons, a lognormal distribution for fracture size may be a better choice in this case.

Stripa 3D fracture model
Based on the optimised parameters described above, the 3D fracture model of the Stripa Mine can be crea- ted using the parameters shown in Table 3.An example is shown in Fig. 11 for one realisation of the model generated for a tunnel of 40 m long with a 4 m x 4 m cross-section.
To assess whether the target statistics are reproduced by the model, 40 independent realisations of the model were generated and sampled on a vertical plane (to represent wall mapping) and a horizontal plane (to represent floor mapping) to derive statistics for fracture trace length and fracture spacing, as shown in Fig. 12.The results are shown in Fig. 13.As can be seen, the target statistics are reproduced reasonably well.For fracture trace length, based on the horizontal (floor) sampling plane, the small size traces are under-represented in the model.This is mainly due to the exceptionally large proportion of smaller size traces on the floor mapping of fracture sets 2 and 3. Overall the match between model and target statistics of trace length is satisfactory, given the compromise used in the parameter optimisation process discussed above.For fracture spacing, the model statistics agree extremely well with the target statistics, suggesting that the fracture densities of fracture sets (with size > 0.5 m) are correctly represented in the model.Note that in Fig. 13, the ± 95% confidence intervals based on the simulations are also plotted.As an additional validation based on the simulations, the censoring statistics and 95% confidence intervals for a vertical sampling plane (wall) are: T0 (no censoring) -84% (81% -88%) and T1 (one end censoring) -16% (12% -18%), compared with the target statistics of T0 -86.6% and T1 -13.4%.The censoring statistics for the horizontal sampling plane (floor) are: T0 (no censoring) -80% (77% -85%) and T1 (one end censoring) -20% (15% -22%), compared with the target statistics of T0 -70.2% and T1 -29.5%.The censoring statistics from the model match very well those of the wall mapping, while they are reasonably close to those of the floor mapping.
To evaluate further the reproduction of target statistics for individual fracture sets, comparisons between model and target statistics are summarized in Figs. 14 -Fig.18.As can be seen from these figures, the target statistics of fracture set 1 are reproduced very well, particularly for the fracture spacing.For fracture set 2, the statistics from wall mapping are modelled reasonally well, but the trace length of the floor mapping is under-represented by the model for smaller values.This is a consequence of the compromise of selecting a l value based only on the statistics of wall mapping, due to the difficulty of finding a suitable value from the fracture traces obtained from the floor mapping, as discussed above.For fracture sets 3 and 4, the match between model and target statistics is in general satisfactory, although the 95% confidence intervals are large due to th e small number of fractures and erratic variations in the derived statistics.For fracture set 5, the model statistics again agree well with the target statistics, particularly for the fracture spacing.
Based on the comparison analysis, it is reasonable to conclude that the simple technique described in this paper is an effective non-parametric method for rock fracture modelling.The target statistics used in this paper are mainly the trace length and fracture spacing distributions, however, the technique is gen-eral and flexible as additional optimisation criteria (e.g., for fracture orientations or fracture connectivity measures) can easily be incoporated into the optimisation process.Note that the work described here is only to demonstrate the framework of the non-parametric method and by no means it can be claimed that the derived fracture model is the "best" model to represent the Stripa data set.As discussed above, the assumption of an exponential distribution for the fracture size may not be suitable in this case.The assumption of a uniform distribution for fracture locations may also be an over-simplification for some fracture sets in this dataset.It is worthwhile exploring other potential models in future work to obtain a better understanding of the fracture model for the Stripa Mine.

Conclusions
This paper describes a non-parametric method to deriving a discrete fracture model based on fracture mappings of multiple sampling planes at different orientations.It differs from traditional approaches in that no   assumption is needed for the distributions of sample statistics, such as the fracture trace length and fracture spacing, and no sampling bias correction is necessary prior to the parameter fitting process.This alleviates the inconsistency problem that sample statistics from different sampling planes for the same fracture set can give different distribution characteristics, as demonstrated in the Stripa dataset used in this studty.
c 2 statistics are used in this approach to quantify the goodness-of-fit of model statistics to the target statistics.To ensure that the comparison is done on the same statistical features, the sampling planes and sampling techniques used for the model should be identical to those used to collect the dataset.The non-parametric approach then searches the parameter space to find the set of parameters that give the best goodness-of-fit of model statistics to the target statistics.To account for the stochatistic nature of the approach, a certain number of Monte Carlo realisations are used to derive the statistics for each set of parameters.Using multiple independent simulations also helps to quantify the uncertainties associated with the model statistics, hence the uncertainties of the derived parameters.
For the Stripa data set, the fracture mappings of the two walls and the floor of the ventilation drift are used as a case study to demonstrate the quality of the described method.Fracture set classification is based on stereographic analysis, which allowed the identification of five sets of fractures including four orientation clusters and one background set with completely random orientation distribution.On the assumption of an exponential distribution for the size of equidimensional 3D fractures, the application of the parameter optimisation based on trace length gives the parameter l values of 2.2 for fracture sets 1 to 4 and 2.8 for fracture set 5. For fracture density, a uniform distribution of fracture locations is assumed and the application of the parameter optimisation based on fracture spacing gives values of the P30 parameter of 0.6 for fracture sets 1 and 3, 1.0 for set 2, 0.8 for set 4 and 2.6 for set 5. The final fracture model for the ventilation drift is constructed on the basis of this set of parameters.The combined model is then sampled and the corresponding statistics are compared with those of the fracture mappings to validate the model.Note that the statistics of all fractures together are not used directly in the parameter optimisation process.The case study demonstrates that the non-parametric approach described here is effective and simple and the target statistics are well captured by the model.
The 3D fracture model derived in this work for the Stripa dataset is based on the assumption that 3D fractures are equidimensional discs with the size distribution following an exponential distribution with only a single parameter l, and fractures of each set are uniformly distributed in three-dimensional space.Thus it cannot be claimed that the model described is the "best" for the Stripa Mine not least because all these assumptions can be challenged.For example, the exponential distribution for fracture size is used in this study because of its simplicity, but it may not be suitable as no significant small size fracture traces were reported in the literature for the Stripa Mine.There is also some evidence that the fractures of some fracture sets may not be equidimensional due to their significantly different signatures of fracture traces on different sampling planes.The Poisson distribution is also used for simplicity but fracture distributions of some fracture sets are clearly not uniform and a non-homogeneous distribution model may be more appropriate.These challenges remain important issues to be addressed in future work.

Figura 8 .
Figura 8. Comparison of target statistics between different fracture sets (horizontal axis is trace length (left) or fracture spacing (right) in metres and vertical axis is cumulative frequency).Figura 8. Comparación de los estadísticos objetivo entre diferentes conjuntos de fracturas (el eje horizontal es la longitud de la traza (izquierda) o espaciado de fracturas (derecha) en metros y el eje vertical es la frecuencia acumulada).

Figure 10 .
Figure 10.c 2 values of the fracture sets at different P30 values (horizontal axis is P30 and vertical axis is c 2 ).Figura 10.Valores c 2 de los conjuntos de fracturas para diferentes valores P30 (el eje horizontal es P30 y el eje vertical es c 2 ).

Figure 12 .
Figure 12.An example of fracture traces of vertical (top) and horizontal (bottom) sampling planes of the DFN model.Figura 12. Un ejemplo de trazas de fracturas de los planos de muestreo vertical (figura superior) y horizontal (figura inferior) del modelo DFN.

Figure 13 .
Figure 13.Comparison of model and data statistics of trace length (top) and fracture spacing (bottom) of all fracture sets (horizontal axis is trace length (top) and fracture spacing (bottom) in metres and vertical axis is relative frequency).Figura 13.Comparación entre el modelo y estadísticos de los datos de longitud de la traza (parte superior) y espaciado de fracturas (parte inferior) de todos los conjuntos de fracturas (el eje horizontal es la longitud de la traza (parte superior) y espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa.

Figure 14 .
Figure 14.Comparison of model and target statistics for fracture set 1 (horizontal axis is trace length (top) and fracture spacing (bottom) in metres and vertical axis is relative frequency).Figura 14. Comparación de los estadísticos del modelo y objetivo para el conjunto 1 de fracturas (el eje horizontal es la longitud de la traza (parte superior) y el espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa).

Figure 15 .
Figure 15.Comparison of model and target statistics for fracture set 2 (horizontal axis is trace length (top) and fracture spacing (bottom) in metres and vertical axis is relative frequency).Figura 15.Comparación de los estadísticos del modelo y objetivo para el conjunto 2 de fracturas (el eje horizontal es la longitud de la traza (parte superior) y el espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa).

Figure 16 .
Figure 16.Comparison of model and target statistics for fracture set 3 (horizontal axis is trace length (top) and fracture spacing (bottom) in metres and vertical axis is relative frequency).Figura 16.Comparación de los estadísticos del modelo y del objetivo para el conjunto 3 de fracturas (el eje horizontal es la longitud de la traza (parte superior) y el espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa).

Figure 17 .
Figure 17.Comparison of model and target statistics for fracture set 4 (horizontal axis is trace length (top) and fracture spacing (bottom) in metres and vertical axis is relative frequency).Figura 17. Comparación de los estadísticos del modelo y del objetivo para el conjunto 4 de fracturas (el eje horizontal es la longitud de la traza (parte superior) y el espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa).

Figure 18 .
Figure 18.Comparison of model and target statistics for fracture set 5 (horizontal axis is trace length (top) and fracture spacing (bottom) in metres and vertical axis is relative frequency).Figura 18. Comparación de los estadísticos del modelo y del objetivo para el conjunto 5 de fracturas (el eje horizontal es la longitud de la traza (parte superior) y el espaciado de fracturas (parte inferior) en metros y el eje vertical es la frecuencia relativa).

Table 1 .
Orientation statistics for the classified fracture sets shown in

Table 3 .
Table3summa-Optimised fracture set parameters based on the assumptions of exponential distributions of fracture sizes and a Poisson distribution for fracture locations.Tabla 3. Parámetros optimizados para los diferentes conjuntos de fracturas basados en las asunciones de distribuciones exponenciales de tamaños de fracturas y distribución de Poisson para las localizaciones de las fracturas.