From outcrop scanlines to discrete fracture networks, an integrative workflow

Understanding fractures and fracture networks is essential for the investigation and use of subsurface reservoirs. The aim is to predict the fractures and the fracture network when there is no direct access to subsurface images available. This article presents a universal workflow to numerically compute a discrete fracture network by combining the 1D scanline survey method, processed with the newly written SkaPy script, together with the multiple point statistic method (MPS). This workflow is applied to a potential geothermal site in Mexico called Acoculco. We use Las Minas outcrops and quarries as surface analogues for the Acoculco reservoir, as Las Minas and Acoculco are both formed by the influence of a plutonic intrusion into the Jurassic–Cretaceous carbonate sequence of the Sierra Madre Oriental in the Trans-Mexican volcanic belt (TMVB). The intrusion is associated with contact metamorphism and metasomatic phenomena, providing the basis for the mining activities at Las Minas. The results obtained using this workflow demonstrate the feasibility of the approach, which presents a solution combining the efficiency of data processing and an interpretation-driven approach to build realistic discrete fracture networks. This workflow can be used in the process of estimating the permeability of a fracture controlled reservoir, with using only scanline surveys data as input. This is essential in the process of evaluating the feasibility to develop an enhanced geothermal system.


Introduction
The use of the subsurface and the exploitation of subsurface resources require prior knowledge of fluid flow through fracture networks, particularly in low permeability rocks. For nuclear waste disposal, for the enhancement of hydrocarbon recovery from a field, or the development of an enhanced geothermal system (EGS), it is fundamental to constrain the fractures and the fracture network. The question we address here is how to predict the fractures and the fracture network of an area when there is no large scale information available, for example from seismic data or aerial images (Unmanned Aerial Vehicle -UAV, or drone). To overcome this potential problem, we developed a workflow to numerically compute a discrete fracture network based on the combination of the scanline survey method performed at single outcrops, processed with the newly written SkaPy B. Lepillier et al.  (2) m in, m out: these stands for ''meter in and out'' as to describe the interval treated in this row; (3) Dir, Surf Dir, Surf Dip : Survey orientation ('Dir' = Direction of the outcrop relative to the North 'N') and slope of the outcrop; (4) Type: stands for the type of structural feature; (5) Str, F Az, F Dip, F Quadrant: ('Str' = Structural measurement direction relative to the North 'N') are used to fill information relative to the structural feature, such as strike, dip, quadrant towards which it dips; (6) F F: is the ''a priori'' structural set; and (7) Q: stands for the number of these features, when a same feature is being repeated within the same interval. reservoir bulk permeability, which includes faults and fractures and the pore network of the rock.
In nature, fractures are organized as networks, from microscopic to regional scale (e.g. Gillespie et al., 1993;Zhang et al., 2016). These networks generally present a regular arrangement (i.e. orthogonal networks, conjugated network; Healy et al. (2015)) which may vary in terms of orientation, size, density and topology  at various scale (from metre scale to reservoir scale). Fracture network geometries vary spatially due to local variation of the stress field and over time. Hence, it is difficult to predict the geometry of these networks. Understanding fracture network starts with the characterization of the fracture geometry. The International Society for Rock Mechanics (ISRM) proposes a method for this fracture characterization (ISRM, 1978). The key elements used to describe fractures are: • the orientation (strike and dip); • the length (also called fracture ''trace'', or ''extent''); • the connectivity (also called ''abutment'' or ''node'' as for example in term of modelling (Peacock et al., 2016)); • the aperture (''mechanical aperture'' when measured at the outcrop); • the filling (cement or clay filling the fracture void); • and finally, the shape, commonly described with a joint roughness coefficient (JRC) (Barton and Choubey, 1977;Tse and Cruden, 1979;Li and Zhang, 2015).
Because the reservoir is not accessible, these fracture networks cannot be measured directly. Therefore, reservoir properties can best be studied at suitable analogue outcrops (Barbier et al., 2012). Analogues are rock sections, made of comparable composition and geometry, cropping out at the surface. Therefore, outcrop analogues can be used to measure the fracture geometries. As explained in Li et al. (2018), an analogue can always have its own local variations from the subsurface reservoir. For this reason, the surface analogue is used as the base case, which can evolve as new information becomes available, to better represent reservoir geometries.
When working on large-scale outcrops, a common method is to use aerial images of wide, well-exposed surfaces, where fracture networks can be explicitly characterized. Unfortunately, there is no suitable outcrop available to be imaged near Acoculco. For that reason, we used another method, based on simple and systematic measurements. The B. Lepillier et al. Fig. 3. Illustration of the method used to build the scanline in SkaPy: Upper part: shows the simplified fractures (as grey straight lines) in the real outcrop (brown background). Red and white dots with numbering from 0 to 8 represent the stations. These stations delimit intervals (here marked in red); Middle part: gives the expression used to calculate fracture positioning from the dataset to computed scanline; Lower part: shows the computed scanline. Small coloured oblique lineaments in between the fractures describe equal spacing.
scanline survey method is a common approach to obtain the statistics required to predict the fracture geometry of the network at the scale of a geothermal field. The method consists of reporting, measuring and describing all fractures visible at the surface of the rock and intersecting with the measuring tape (ISRM, 1978;Lavenu et al., 2014). The dataset, created from the scanline survey, can then be used for computing a global discrete fracture network (DFN).
A Discrete Fracture Network (DFN) model is a mathematical representation of fracture distribution in space and fracture characteristics. As explained in Lei et al. (2017), the DFN concept is broad, since the geometry of the DFN depends on the need for which it is created, or on the way the fracture characteristics are obtained. For example, a DFN can be generated by digitizing fractures on large scale imagery or it can also be obtained from samples at a laboratory scale. In the context of this article, we compute the DFNs numerically from scanline datasets. We use the scanline datasets as input to extrapolate to a much larger geographical domain. We assume that the fracture parameters measured at the outcrop are representative of the geometry of the network in its near neighbourhood. The extrapolation far from this point is unknown and has been taken into consideration by means of gradual transitions from one sampling area to the other. Common approaches as used in standard software tools would achieve this using a stochastic method, which honours the scanline input data and returns statistically accurate results. In this workflow, we extrapolate the scanline surveys with the multiple point statistics (MPS) method, which in addition to being statistically accurate preserves the geological patterns and therefore preserves the geometry of the fractures.
This article presents a new workflow using the specially developed Python TM2 based SkaPy script, which computes statistics on multiple 2 Python Software Foundation (2019). scanline surveys, classifies the fractures, and represents them in a georeferenced system. This representation is further analysed to provide the input for the calculation of the DFN by the MPS method proposed by Bruna et al. (2018Bruna et al. ( , 2019. In this article, Bruna et al. (2019) analyses applicability of the MPS method to reproduce a known fracture network (i.e. manually interpreted from an outcrop), preserving the geometrical characteristics of the fractures. This workflow combines the efficiency of data processing and an interpretation-driven approach. The script is freely available on GitHub TM .
We apply the present workflow to the Acoculco geothermal field, within the framework of the GEMex project (e.g. Jolie et al., 2018), where two exploratory wells reached high temperatures (300 • C) at reasonable depth (2 km). Unfortunately, the wells did not access any productive geothermal layers. For this reason, Acoculco could eventually be developed as an EGS. Thus, the objective of this paper is to provide a characterization of the existing fracture network prior to stimulation.

Existing methods for the scanline survey
The linear scanline survey method has been widely described (ISRM, 1978;Lavenu et al., 2014;Watkins et al., 2015). The principle of this method is based on laying a measuring tape along the outcrop and reporting every fracture crossing this tape, describing the fracture characteristics such as strike, dip, length, connectivity, aperture, filling and shape. Even though the method is known for potentially biasing the fracture representativity (Terzaghi, 1965), it is still broadly used. Several studies have discussed the fracture length and the fracture representativity in the scanline survey (Priest and Hudson, 1981). As more recently explained in Zeeb et al. (2013)  length. Therefore, short fractures are underrepresented'' (p. 13). This under-representation impacts the statistical distribution of the fracture lengths. The measured fracture lengths are also controlled by the outcrop dimensions. Nonetheless, the method remains a good tool thanks to its simplicity and reliability to give exhaustive characterization of the fractures. Some variations of this linear scanline survey exist as for example the circular window surveys proposed by Mauldon (1998), Mauldon et al. (2001) and a mixed method of transect linear scanlines and circular windows Watkins et al. (2015).

Existing scripts for scanline or fracture data processing
Markovaara-Koivisto and Laine (2012) published a MATLAB ®3 script for scanline data processing. The script automates the dataset computation: it classifies the fractures into sets, provides statistics on the fractures, and builds a 2D or 3D visualization of the scanline. The advantage of their MATLAB ® script is that it handles the shape of the fractures, comparable to the joint roughness coefficient (JRC), named ''undulation'' of the discontinuities. FraNEP is a script developed in Visual Basic TM , written by Zeeb et al. (2013). The FraNEP script automatically analyses the statistical properties of 2D fracture networks by applying the linear and circular scanlines methods. The script includes length and orientation corrections in the calculations and also analyses fractures as a network, calculating their spacing, intensity and length distributions. The script includes the fractures classification into sets and computes rose diagrams of the fractures. More recently, Healy et al. (2017) published a MATLAB ® script under the name FracPaQ. This script is intended to be a ''software-based toolbox'' to quantify fracture patterns in 2D. Similarly to FraNEP, the FracPaQ toolbox is mainly implemented to process photography or interpreted imagery of fracture networks and to analyse the statistical distributions of the fracture properties. Among other results, the code estimates the bulk 3 The MathWorks Inc. (2019). permeability of the system studied. The strength of FracPaQ is that it can process a wide range of spatial scales.
None of these approaches combine a multiple scanline datasets analysis, classification of the fractures, and geographic representation. To transfer this information into a DFN, we needed to develop our own method. SkaPy is designed to fill this gap.
One advantage of scripting is to automate a workflow, which increases efficiency and reduces personal bias. This means, scripting also contributes to decreasing the risk of error. Another advantage is the possibility to analyse the statistical distributions of the fractures. Even more important is the possibility to run these statistics on multiple scanline surveys together. However, to compute any statistics on the fractures, it is fundamental to first separate the fractures into sets.
SkaPy has been written in the continuation of the Markovaara-Koivisto and Laine (2012) script, which processes scanline datasets. SkaPy is a flexible tool that handles surveys with constant or varying orientations. This allows surveys along curved outcrop surfaces. To compute the statistical analysis on fracture characteristics, SkaPy analyses not only one scanline dataset but processes multiple datasets. Therefore, it becomes easy to constrain the fracture characteristics over multiple outcrops. Then, as in the Markovaara-Koivisto and Laine (2012) script, SkaPy generates statistical plots, such as histograms, box-plots and stereonets. Because visualization is also important, SkaPy plots the scanline(s) in a geo-referenced system. A geo-referenced representation of the discretized fractures is the first step towards extrapolation to a DFN.

Existing methods for DFN computation
Many tools already exist and are still being developed to improve DFN generation. Commonly used programs and algorithms for this purpose include the Move TM4 software from Midland Valley, and     FracMAN TM5 from Golder Associates, which was implemented in Petrel TM6 (Schlumberger). FracMAN TM is very effective at assessing the influence of fractures on bulk permeability (Chesnaux et al., 2009). The two most common methods to express the role of fractures on the reservoir flow are: (i) conversion of the DFN into an equivalent 5 Golder Associates (2019). 6 Schlumberger Limited (2019). porous medium (Dershowitz et al., 2004); (ii) or more generally, to upscale the DFN into a continuum of fluid flow model representative of the fractured reservoir (Surrette, 2006). For creating DFNs, most of the methods, such as FracMAN TM or Move TM , use similar approaches based on statistics to populate a simulation domain. The results of these models are statistically correct but often do not represent the geology correctly. The multiple point statistic approach is an emerging method that, in addition to being statistically accurate, identifies the fracture network patterns and preserves them throughout the extrapolation.

Method: Workflow description
The workflow developed in this article can be summarized into four steps ( Fig. 1): (1) the data collection using the linear scanline survey. The data, from one or several surveys, is gathered in one table-like file format (such as ''.csv''); (2) the file is loaded in the Python TM script SkaPy which processes the dataset; (3) using the output from SkaPy, the users manually create, so called Training Images (TIs), an extrapolation of the fractures and their distribution to the outcrop scale and a probabilistic map; (4)The MPS method computes the DFN.

Using SkaPy
SkaPy is made of two scripts. The first one classifies, analyses and plots the scanline; the second one assigns the geo-referenced coordinates. To run properly, SkaPy uses common Python TM open-source libraries such as Numpy, Matplotlib, Pandas and mplstereonet. The input file is loaded and converted into a Pandas DataFrame, which is a Python TM library implemented to manipulate numerical tables. For the analysis of the scanlines, the input data must contain at least the following information: the name of the outcrop; the name of the survey; the lithology; the position along the scanline survey; the orientation of the survey and the slope of the outcrop; the type of discontinuity (stylolite, fracture, fault, dyke); the structural measurements (strike, dip); the ''a priori'' fracture family, even though this is re-calculated during the script compilation; and the quantity ''Q'' that this fracture type has been observed in the interval; the aperture of the fracture (and when possible the filling material) and the length of the fracture trace. To be able to plot the scanline in a georeferenced system, the geographic coordinates of the starting point of the scanline survey are necessary (Fig. 2).
The first SkaPy script is split in seven parts. In SkaPy part−1, the script plots the statistical distribution and histogram of the length and aperture measurements for the whole dataset. This gives a global overview of the data to identify eventual correlations between the two parameters, the fracture trace length and the fracture aperture. In SkaPy part−2, the users need to define the values to be used for the classification of the fractures into sets. The script returns a stereoplot with all the structural measurements, and then one stereoplot per set defined. Proceeding on from this classification, the script computes boxplots, showing the distributions of the fracture heights and fracture apertures per set. The box-plot format emphasizes which values are representative of the core of the distribution as opposed to the ones behaving as outliers. (i) In SkaPy part−3, a first batch of plots is generated, gathering the data per outcrop. (ii) In SkaPy part−4 a second batch of plots is generated, gathering the data per lithology. In SkaPy part−5, the script builds and plots the scanlines. Fig. 3 shows the method used to script the construction of the scanline. For clarity, the fields ''m in'' and ''m out'' of the table are here referred to as ''stations''. The offsets calculated between these stations would be called the ''intervals''. In a first loop, SkaPy builds the trace of the scanline survey: it calculates the length of an interval by determining the positions of the stations, and compiles these intervals following the orientation given in ''Surf Dir''. Then, in a second loop, SkaPy calculates the quantity of fractures Q belonging to each of the fracture sets within this interval. The algorithm then splits the intervals into ( + 1) sections in order to equally distribute the fractures, and to assign a specific position to each of the fractures. X and Y coordinates are calculated from the origin station of the scanline.
In SkaPy part−6, the script is used for the extrapolation of the scanlines into training images (TIs). In this part, instead of using the fracture length and strike values as in SkaPy part−5, the script plots the average strike value representative for each set, and extends the lengths of the fractures such that they fill the entire plot. The objective of this step is to offer guidance to the users for the following step of the workflow, which consists of creating the training images (TIs). In SkaPy part−7, the script computes the stereonets and rose diagrams for each outcrop. The second script, SkaPy−UTM, assigns the UTM coordinates to the scanline, including all the fractures, by modifying the origin of the reference system.

Method: DFN creation using Multiple point statistics
Multiple point statistics (MPS) is used to extrapolate geological patterns to larger scales. Classical geostatistical methods use variograms; the MPS method uses the training images (TIs) as input data for the simulation. The concept is based on analysing the pixel content of the TI to identify the spatial patterns to be reproduced into a larger domain (Tahmasebi, 2018). The TI is considered as a grid of pixels which contains geological patterns. For that reason, the TI needs to include the possible range and shape of the geological entities to be modelled . In this study, we used the Direct sampling (MPS) code as presented in Straubhaar et al. (2011).

Creation of the training image
The TIs of the fracture network, created from the scanline generated by SkaPy (Fig. 4), are then extrapolated into a larger model, creating the DFN. Scanline surveys provide datasets of very localized data, usually from tens to hundreds of metres. In comparison, the DFN is often used to represent fracture geometries over a much larger area, such as kilometre scale reservoirs. Building the TIs is the opportunity B. Lepillier et al.

Lithology
EAC 1 well-marker depth EAC 2 well-marker depth SkaPy, the georeferences are derived from the scanline coordinates.

Training images and probability map
The TI is the only input for fracture geometries in the MPS method. For this reason, all ranges of dimensions, orientations of the fracture patterns need to be represented in the TI. The TI is saved as an image format, so that it can be analysed as a grid where every pixel has a colour value coding the structural entity. A reference table is given in Table 1 as an example. The algorithm of the MPS analyses the value of every pixel and pixel associations to identify the pattern composing the TI.
The MPS then builds one DFN based on two TIs and a probability map (Fig. 5). The probability map is a representation of the destination B. Lepillier et al. Fig. 13. SkaPy output part−1: The dataset at glance, from left to right: exceedance frequency of the fracture heights, then exceedance frequency of the fracture apertures and finally plotting fracture apertures against fracture heights; fractures are in this last plot sorted by using fracture numerical index from the database. As the database consists of a succession of scanline surveys, fractures are plotted in the following order: (from start to end of surveys) Boquillas, Eldorado, Pueblo Nuevo, Tatatila, Rinconada, San Antonio Tenextepec. domain where to populate the extrapolated fracture networks. This map defines the areas where each TI should be used to populate the DFN in the simulated domain. The probability map could represent geological entities, such as faults or fluvial channels or different lithologies. In this study, however, we linearly interpolate between the two TIs for simplicity.

Acoculco within the regional geological context
The method we developed creates discrete fracture networks (DFNs), where the main input are scanline datasets measured at outcrops near the Acoculco test site.
When reporting a scanline in the field, measurements are done objectively. Later, when processing the data, the users want to set these measurements in context. This is why it is essential to know the regional structural context to classify the scanline dataset.
Acoculco is located in the Trans-Mexican volcanic belt (TMVB), about 100 km North-East of Mexico City. The TMVB is a mountain range resulting from the subduction of the Avalon, Rivera, Orozco and Cocos plates beneath the North American plate (Ferrari et al., 2012;Manea et al., 2013) (Fig. 6). The TMVB is composed of a thick series of Pliocene to Quaternary volcanic deposits forming a high volcanic plateau. This plateau overlies a fold and thrust belt resulting from the Sierra Madre Oriental, a Laramide deformation of the Jurassic to Cretaceous deposition of the carbonate sequence (Carrasco Núñez et al., 2017;Norini et al., 2015;López-Hernández et al., 2009;Campos-Enriquez and Garduño-Monroy, 1987).
Two geothermal exploration wells, EAC1 & EAC2, were drilled in Acoculco, about 500 m apart, in 1995 and 2008. Since the wells are really close to each other and even though the well-markers position suggest a fault could cause an offset between the two wells, their lithological sequence is very similar and can be summarized, as described by (Viggiano-Guerra et al., 2011;López-Hernández et al., 2009), in Table 2: Fig. 7 is a representation of the geological section in Acoculco adapted from the model of López-Hernández et al. (2009), based on gravity, magnetic and well data, which has been further developed by Lorenzo Pulido et al. (2010) and Canet et al. (2015). The carbonate section consists of dolomitic limestones, limestones metamorphised into marbles or metasomatised into skarns due to the intrusion of a granodioritic pluton.

Las Minas analogue
Las Minas is a village located down in a valley, about 120 km East of Acoculco. Las Minas valley is considered here as an excellent analogue of the Acoculco geothermal system. As seen on the geological map (Fig. 8), Acoculco and Las Minas are made of rocks that were deposited B. Lepillier et al. in the same geological environment. The exposed rock sequence in the Las Minas analogue is comparable to the log of Acoculco wellbores, from the recent volcanic cover, overlaying the carbonate sequence made of the Early Cretaceous limestones called the Tamaulipas Formation, all the way to the granodioritic pluton, which is exposed deep in the valley, where the river Rio Las Minas runs. Hundreds of metres away from the river bed, two outcrops, named Boquillas and Eldorado, expose skarns; higher on the flanks of the valley, two marble outcrops are found: Pueblo Nuevo and Tatatila. Further away is first one outcrop made of marblized dolomitic limestones called Rinconada, and much further, beyond the influence of the intrusion, San Antonio Tenextepec is made of unaltered dolomitic limestones. These outcrops represent a proximal to distal context for the granodioritic intrusion into the carbonates. Fig. 9 is a simplified representation of the geological context.

Structural context
Four major structural trends are identified in the region from Acoculco to Las Minas. The Acoculco Caldera is located at the intersection of two regional fault systems, one trending NE and the other NW (López-Hernández et al., 2009). The main structures in the region are inherited from: (i) the Laramide Orogeny, Late Cretaceous with a NE -SW compression, inducing NW -SE thrusts and folds affecting the lithological sequence from the metamorphic basement up to the Mesozoic sedimentary rocks; (ii) the Eocene-Pliocene extensional and transtensional phases that produced N -S to NE striking faults (Carrasco Núñez et al., 2017;Norini et al., 2015). A micro-structural study identified two major structural directions: (i) the first direction is from N140 to N170; (ii) these structures are cut by a second trend oriented from N40 to N70; this confirms the regional trends and identifies the direction N40 as a maximum regional stress (Campos-Enriquez and Garduño-Monroy, 1987).
In the Los Humeros Volcanic Complex (a ten of kilometres westwards from Las Minas), the main structures driving the ascending hot fluids have been identified as the NNW -SSE east-dipping faults (Norini et al., 2015). Additionally, Carrasco Núñez et al. (2017) report that producing geothermal wells are located along the main NNW -SSE active faults or near the N-S striking fault splays NW -SE and NE -SWtrending regional faults systems appear to have been permeable fluid flow pathways (López-Hernández et al., 2009).

Case Study: Applying the scanline survey
In the Acoculco field, we want to predict the distribution of the fractures in the reservoir formation which corresponds to the carbonate sequence. As described in Section 4.1.2, the carbonate sequence includes three lithologies: the limestones, marbles and skarns. For that reason, the linear scanline survey method was recorded in two outcrops for each of these rock types, for a total of six outcrops (Table 3).
The trace of the scanline survey follows the direction of the outcrop walls. As explained in Section 3.1, each interval is delimited by two stations in between which the wall follows a constant direction (Fig. 10). Within these intervals, fractures are counted individually or per group of fractures having the same strike and dip values; the number of fractures is then reported in column ''Q'' of the dataset (Fig. 2). The value reported for the aperture is the average value of the aperture measured for a fracture, or a set of fracture. The length and height of the fractures are measured along the outcrop. Because of that, if a fracture trace is longer than the outcrop, its measurement is limited to the length of the trace visible within the outcrop. The fracture trace length measurements closely depend on the morphology of the outcrop. The height and slope of the outcrop walls change from the start to the end of the survey. For that reason, we applied a correction on our fracture trace length measurement, by converting it into a vertical height and a horizontal length, based on the slope of the outcrop, using a simple trigonometric rule (see Fig. 11). The measurements are well constrained up to 5 m (±1 × 10 −3 m), fairly well constrained up to 10 m and eventually less over 20 m (±5 m), depending on the regularity of the outcrop morphology and its accessibility for the measurements. Fig. 12 shows a sample of the Las Minas outcrops.

Case Study: Running SkaPy
The dataset, containing the scanlines from the six outcrops, is composed of about 600 rows. The total length of scanlines is 257 m, which includes 2090 fractures (Table 4).
According to this dataset, most of the fractures in the skarn outcrops, Boquillas and Eldorado are tight to closed (in this study, fracture aperture classification results from qualitative and quantitative description. As indication, open fracture referrers to aperture > 0.1 cm, partly open: aperture of the fracture varies from 0.05 to 0.1 cm, tight to partly open: aperture of the fracture varies but is always < 0.05 cm, closed: the fracture never has any aperture). In the marble outcrops, the quantity of fractures is much lower, but most of them are open to widely open. Regarding the limestone outcrops, Rinconada is located on a regional fault and San Antonio Tenextepec is on a thrust fault. Both of them are highly fractured. As the scanline in San Antonio Tenextepec is the longest, a higher range of variability within the dataset is expected.
When running SkaPy, the first figure shows the distribution of fracture heights and fracture apertures for the whole dataset. In this case there is no direct correlation between the size of the fracture and the width of the aperture (Fig. 13).
In SkaPy, part 2, we first classified the fractures in 4 sets, based on the regional structures described in the literature presented in Section 4.1.3. These 4 sets are: N50, N100, N135, N170. Since there is a large amount of fractures, we decided to apply a more detailed classification, adding 3 more sets, based on the dip angle of the fractures (Fig. 14).
As the objective is to characterize the fractures, we want to calculate the representative values, per lithology, of the fractures height and aperture. In Fig. 15, the left column shows the box-plots and distributions of fracture heights and apertures measured in the outcrops, and the right column shows the box-plots and distributions recalculated for all the examples of a given rock type. Then, SkaPy builds the scanlines, and SkaPy-UTM plots them in the UTM system (Fig. 16).

Case Study: Creating the Training Images (TIs) and running MPS
Following the workflow chart (Fig. 1), the next step is the creation of the TIs. Fig. 4 shows the method used to create the TI of the Boquillas outcrop, one of the two skarn outcrops. In this example, the scanline survey is made of two orthogonal sections separated by 20 m: the main one is about 25 m and the second one 6 m long. We used the output from SkaPy part−6, where the extended fractures serve for guidance. This tool emphasizes the fracture patterns. In this case, the scanline shows that the blue and red sets, representing sets F3 and F4 respectively, are organized as clusters of fractures, while the sets yellow F7 and green F2, are much more equally distributed. This extrapolation works as the process of up-scaling: the scanline is a very high resolution dataset (decimetric), extrapolated to a DFN of a much larger scale, as in this case (600 x 600)m 2 . Therefore, not all fractures can be populated, and the TI serves as a transition tool representative of the scanline survey at a large scale.
Following this method, Fig. 17 shows the TIs obtained for the three lithologies: the skarns, the marbles and the limestones. These training images, together with a probabilistic map, as presented in Fig. 5, are used as inputs for the computation of the DFNs, using the multiple point statistics method.
In this case study, the objective is to obtain a DFN model representing the transition of the fracture distributions, from one outcrop to the second, within the same rock type. For this reason, we used a linear extrapolation from one TI to the other. Therefore, the probability map consists of a domain of (600 x 600) m 2 , divided into three equivalent sub-domains of (200 × 600) m 2 . The first sub-domain is populated using one TI, the second with the other TI and the third, in the middle is populated with the calculated interpolation.

Case Study: created DFNs
The fracture patterns are very well preserved during the MPS extrapolation (Fig. 18). For example, the marble rock type stays, after extrapolation, consistent to its sparse and homogeneous fracture distribution, while the extrapolation of the fracture distribution for the limestones well represents the influence of regional faults, which causes the higher fracture frequency. Another element to control the fracture patterns extrapolation is the fracture lengths. The 3 created DFNs show longer fractures on their right side, corresponding to the influence of the second TIs of these three lithologies.

Discussion & way forward
In the process of generating DFNs, several issues affect the quality of the resulting model. When DFN is created following the method described here, one has to be aware that the error or biases induced by the scanline survey itself, as for example the fracture trace length or fracture representativity, are not overcome. Scanline data collection also present a limitation regarding fracture relationships, as in this study, we do not report information regarding fracture abutments. A first clear solution, would be to complete the scanline dataset with interpretation of virtual outcrop images from drone imagery as developed by NORCE Research with Lime TM software (Buckley et al., 2019). This would mitigate the issue of fractures representativity and provide data on fractures relationships. Eventually, this could be coupled with the Facets plugin developed by Dewez et al. (2016) and implemented in the CloudCompare TM platform. 7 This Facets plug-in is an algorithm that detects automatically the orientations of the outcrop surfaces. Assuming these surfaces are the result of erosion working along weaknesses of the rock, the algorithm could be used to estimate the most frequent fracture directions over an outcrop.
Another solution would be to integrate the circular scanlines surveys at regular intervals along the linear one, as proposed by Watkins et al. (2015), to obtain a more complete dataset and reduce these biases, providing data on fracture chronologies and topology relationships.
Another limitation is the accessibility of the reservoir itself. In this study we used Las Minas as analogue system of the Acoculco geothermal field. This gives us a good prediction of the fracture network at the surface. By essence, there might be differences regarding past and present tectonic stresses between Acoculco subsurface reservoir and its surface analogue of Las Minas. To get a more accurate model, it would be possible to correct Las Minas fracture models by calibrating to the Acoculco well core samples. In addition, the present day local stress field in Acoculco at depth may differ from the one at the surface of Las Minas. This may be compensated by measuring the current stress field downhole, and applying a correction on the fracture aperture and orientation. More generally, using well log data, such as Formation 7 CloudCompare (2019). Imagery (FMI), would significantly refine the prediction of the fracture distribution in the subsurface. Regarding the approach itself, three elements related to SkaPy have to be considered: Markovaara-Koivisto and Laine (2012) in their MATLAB ® script, include an automatic classification of the fractures into sets. SkaPy does not automatize this part; we ask the users to define their own sets. We prefer to leave this classification to the users rather than to statistics. This could be debated as it makes room for personal bias. The choice of leaving this is based on the assumption that the users know which fracture sets to expect, and reinforces the will to keep control on the geological meaning of this processing. The users can control the quality of their classification by checking the plotted stereonets (Fig. 12). Nevertheless, the automated classification tool could be implemented for guidance.
In the SkaPy scripts, as described in this manuscript, fracture geometry is not taken into account. We see two possible ways to improve this limitation: the first would be to integrate the fracture abutments; and the second is the classification of fracture shapes proposed by ISRM (1978).
A starting point to define the fracture abutment could be made by defining a chronology of the fractures determining the 'priority' on fracture intersections.
Regarding the fracture typology, the populating of the TIs, could follow a similar classification, as initiated by Markovaara-Koivisto and Laine (2012) describing the fracture tortuosity.
In SkaPy part−5, at the scanline construction, the script induces a small error of fracture positioning, which is described in the Fig. 3. As explained in 3.1, the scanline is built using ''stations'' and ''intervals''. In these intervals, fractures from the same set are distributed at equal distances over the interval. As a consequence of that method, a single fracture will be placed in the middle of the interval instead of its exact position in reality. This error is kept small by the choice of interval length. The larger the interval, the bigger the error. This workflow has been developed to analyse the feasibility of establishing an enhanced geothermal system in the Acoculco field.  For that reason, we need to predict the fractures and the fracture network structuring the reservoir formations. The aim is to be able to calculate the role of the fractures for the fluid flow circulation. The DFN provides the basis for simulating the fluid circulation using a doublet of two wells, one injector and one producer, connected by fractures. Preliminary results are presented in Fig. 19 and the full study is available in Lepillier et al. (2019).

Conclusion
In the context of developing an enhanced geothermal system (EGS), it is fundamental to evaluate the fracture system present in the subsurface. Very often no detailed information for such an evaluation is available. In the study presented here, we propose a workflow to predict the fractures and the fracture network from outcrops. The four steps from data collection in suitable analogue systems at the surface to generating a Discrete Fracture Network provide results that go well beyond a statistical representation of the fracture networks observed at the surface. By preserving the real fracture distribution as well as measured fracture characteristics such as apertures, derived lengths and orientations, the DFNs computed this way can be used to predict reservoir properties, such as bulk permeability and provide the basis for reservoir enhancement procedures for a potential EGS development.
While the method cannot replace information retrieved from direct access to the subsurface by downhole images and measurements, it does provide, to our knowledge, the currently most realistic extrapolation of a surface fracture model from outcrop to field scale.

Data availability
The dataset used and the SkaPy script presented in this study are available in the repository: https://github.com/BatLep/SkaPy.git.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.