FracPaQ: A MATLABTM toolbox for the quantification of fracture patterns

The patterns of fractures in deformed rocks are rarely uniform or random. Fracture orientations, sizes, and spatial distributions often exhibit some kind of order. In detail, relationships may exist among the different fracture attributes, e.g. small fractures dominated by one orientation, larger fractures by another. These relationships are important because the mechanical (e.g. strength, anisotropy) and transport (e.g. fluids, heat) properties of rock depend on these fracture attributes and patterns. This paper describes FracPaQ, a new open source, cross-platform toolbox to quantify fracture patterns, including distributions in fracture attributes and their spatial variation. Software has been developed to quantify fracture patterns from 2-D digital images, such as thin section micrographs, geological maps, outcrop or aerial photographs or satellite images. The toolbox comprises a suite of MATLABTM scripts based on previously published quantitative methods for the analysis of fracture attributes: orientations, lengths, intensity, density and connectivity. An estimate of permeability in 2-D is made using a parallel plate model. The software provides an objective and consistent methodology for quantifying fracture patterns and their variations in 2-D across a wide range of length scales, rock types and tectonic settings. The implemented methods presented are inherently scale independent, and a key task where applicable is analysing and integrating quantitative fracture pattern data from micro-to macro-scales. The toolbox was developed in MATLABTM and the source code is publicly available on GitHubTM and the MathworksTM FileExchange. The code runs on any computer with MATLAB installed, including PCs with Microsoft Windows, Apple Macs with Mac OS X, and machines running different flavours of Linux. The application, source code and sample input files are available in open repositories in the hope that other developers and researchers will optimise and extend the functionality for the benefit of the wider community. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Background
Fractures and their patterns exert a fundamental influence on the mechanical and transport properties of rocks. Fractures rarely occur in isolation, and their patterns are often highly complex. Length scales can vary from sub-micrometre to many kilometres (Bonnet et al., 2001). Fracture orientation distributions also vary widely, and their spatial densities are often heterogeneous (Gillespie et al., 1993). Understanding fracture patterns e the geometrical attributes of the constituent fractures and of their total ensemble e is important in many sub-disciplines of earth sciences: structural geology and tectonics, impact geology, rock physics and geophysics, hydrogeology, energy and storage of hazardous products (da Silva et al., 2004;Gurov and Koeberl, 2004;Rao et al., 2000;Su et al., 2001;Watkins et al., 2015b). The quantification of an observed fracture pattern is a necessary precursor to understanding the physics underlying its formation, and to making robust predictions about its extent and scaling in the subsurface, which ultimately determine the transport properties of the network (i.e. permeability).
Scope clearly exists for practitioners, including both academia and industry, to use a common set of tools to measure and quantify patterns of fractures in rocks (e.g. Hardebol and Bertotti, 2013;Zeeb et al., 2013a). For example, DigiFract is a software package written in Python designed to work directly with fracture data digitised from outcrops, and is based around a geographical information system (GIS) core (Hardebol and Bertotti, 2013). FraNEP is an Excel™ based tool written in Microsoft Visual Basic™ for analysing fracture networks using a range of spatial sampling methods (Zeeb et al., 2013a).
Other programs for processing data acquired by laser scanner and borehole televiewer measurements also provide some basic tools for the analysis of fracture lengths and strike (e.g. Masoud and Koike, 2011;FracMan7, 2012). Markovaara-Koivisto and Laine (2012) provide a MATLAB script for the analysis and visualization of scan line data. The packages LINDENS (Casas et al., 2000) and SAL (Ekneligoda and Henkel, 2010) use the coordinates of fracture end points, for example from GIS analysis (e.g. Holland et al., 2009a), as inputs. Both programs analyze fracture length and orientation using frequency histograms and rose diagrams. The first one provides information about fracture density, whereas the second one provides additional properties such as fracture spacing and unidirectional frequency. The software package FracSim3D (Xu and Dowd, 2010) is a fracture network generator, but incorporates scan line, window and planar methods to sample fracture network characteristics. FracSim3D also offers statistical tools including histogram analysis, probability plots, rose diagrams and hemispherical projections (Xu and Dowd, 2010). These examples of existing packages and tools illustrate that such programs are often developed for a specific study or purpose. To our knowledge, no one open-access software provides for the comprehensive and complete analysis of rock fractures in a readily available and easily extensible programming language (e.g. MATLAB).
This paper describes a new collection, or toolbox, of MATLAB™ programs designed to quantify fracture patterns in two dimensions (2D) from digital data (Fig. 1). The input can either be a binary image file of fracture traces (e.g. a. JPG/.JPEG or. TIF/.TIFF image) or an ASCII tab-delimited text file of (x,y) coordinates that mark the nodes of each fracture trace. The outputs from the toolbox include quantitative estimates for the attributes of individual fractures (e.g., their lengths and orientations) and for the attributes of the whole pattern (e.g., connectivity, permeability). The toolbox is currently limited to 2D fracture patterns and the input area is assumed to have no topographical variation. FracPaQ is the first open-source, cross-platform, freely available toolbox based on published methods. We believe that FracPaQ can provide a foundation for a shared set of open and objective methods with the advantages of readily reproducible results, and customisation and extension through the MATLAB scripting language (Fig. 2). This paper is organised as follows: background information about the rationale and issues of quantifying 2D patterns is provided in Section 2. Section 3 introduces the example datasets used throughout the paper. Section 4 provides a description of the program structure, inputs, processing and outputs, with examples. Section 5 discusses the limitations and scope for further work, and is followed by a short summary in Section 6.

Rationale
The fundamental premise of this paper, and rationale behind the development of FracPaQ, is that a fracture pattern is much more than a set of orientations shown on a stereogram. Orientation distributions are a key component of any fracture pattern, but length distribution, spatial density and connectivity are equally important to defining the fabric of fractured rock, and the "full characterisation of fracture patterns will therefore require independent analysis of each of these attributes" (Gillespie et al., 1993).
Quantitative characterisations of fracture patterns can also be used to estimate important properties such as permeability, rock mass integrity and strength. In addition, it can be important to quantify the statistical distribution of size attributes (e.g., length, aperture) through scaling laws. FracPaQ can be used to quantify fracture attributes and estimate properties using the same input data.
It is well known that 2D fracture datasets have limitations, such as truncation, censoring, and 'cut effects' (e.g., Zeeb et al., 2013b and references therein). Truncation refers to the under-representation of smaller fractures due to the resolution limits of the sampling method (Bonnet et al., 2001). In the case of using digitised fracture trace maps, this limit is effectively the size of a single pixel in the original image. Censoring refers to the effect that arises when the sample window (e.g., and outcrop) is smaller than the trace length of the longer fractures, so that those lengths are lost from the analysis. 'Cut effects' refers to the potential geometric bias introduced by the intersection of variably oriented 3D fracture planes with the 2D plane of analysis.
We have accepted the technical challenges outlined by Bonnet et al. (2001, their Section 8.6 "Future Research"), by developing FracPaQ as a tool to help: find objective ways to analyze simultaneously for orientation and spatial position; take account of the subtle spatial relationships between fractures in a population; quantify the spatial correlation of objects with finite length; and pose the scaling hypothesis to be tested more objectively than is commonly the case.
Fracture patterns in rocks are 3D phenomena, and yet generation of fully 3D fracture pattern datasets is quite difficult. Highresolution reflection seismic volumes, LiDaR datasets from large outcrops and X-ray CT scans of core plugs can provide some insights into three-dimensional pattern geometry as a function of resolution, penetration depths and quality of contrasts in material properties. Consequently, these imaging methods have limitations for constructing a truly 3D dataset. LiDaR datasets have been labelled '2.5D', as they do not provide access to the interior of the scanned outcrop surface. Subtle complications in CT datasets can produce images with poorly resolved geometries, e.g. fracture boundaries may be weakly defined due to artefacts known as the 'partial volume effect', which if not corrected can lead to erroneous quantification of feature dimensions (Ketcham and Carlson, 2001). Fully 3D fracture network datasets are still relatively rare, can be expensive to acquire and are invariably incomplete. In contrast, 2D datasets of fracture patterns in geology and geophysics are almost ubiquitous from satellite and aerial photographs, to outcrop photographs and maps, to thin section photomicrographs (either via optical or electron microscopy). In this paper, we have generated fracture trace maps from images of fractured granite, limestone and sandstone, spanning a total scale range of greater than 100 m to less than 1 mm. FracPaQ can exploit and quantify the abundance of 2D fracture data, with a future aim of the quantifying 3D patterns from multiple 2D datasets.

Program description
FracPaQ is designed to generate quantitative fracture pattern data with user control over the output. The application is installed by copying all of the files from the GitHub (http://davehealyaberdeen.github.io/FracPaQ/) or MathWorks FileExchange repositories into a folder on the user's computer. After starting MATLAB™, the working directory should be set to the folder containing the source code and sample files (Fig. 3). The application is started by typing 'guiFracPaQ2D' in the command window of a MATLAB™ session. The main processes in the flow chart are mapped to similarly named buttons in the graphical user interface (GUI), e.g. Browse …, Preview and Run. The decision to implement FracPaQ in MATLAB™, and not another coding language such as Python or Cþþ, was based on a number of factors, including general availability, and familiarity with MATLAB™ among structural geologists, a wide range of useful built-in functions in the various add-on Toolboxes (e.g. for Image Processing, and Statistics), and relative ease of use.

Input file formats
FracPaQ currently accepts two main types of input file data: 1) tab-delimited ASCII text files of fracture trace nodes ("node file") and 2) graphical image files of fracture traces ("image file"). Supplying a node file of specific (x,y) coordinate pairs of every node along every fracture trace is the most robust way of entering data into FracPaQ. The user is free to generate this input file in their preferred software. For the examples used in this paper, we have prepared node files of fracture trace maps in a 3-stage pre-process: These examples demonstrate the ranges in length scale that can be addressed using the FracPaQ toolbox, from remote sensed maps of 10 2 m to micrographs at less than 10 À3 m (<1 mm).
1. Importing a digital photograph into Adobe Illustrator™ and then manually tracing fractures with the pen tool, which generates 'line' or 'polyline' elements, onto a new layer; 2. Saving the layer with the fracture traces (i.e. without the underlying photograph) as a scalable vector graphics file (.SVG), so that the traces are saved as 'line' and 'polyline' tags; 3. Using a custom C-shell script (svg2fracpaq.csh; included with the FracPaQ source code) to extract the (x,y) fracture-trace coordinates from this. SVG file and write them into a tab-delimited text file.
Our use of a custom C-shell script currently restricts this particular method of pre-processing to those operating FracPaQ on Unix or Mac OS X systems. Unix shell script support may arrive in newer versions of Microsoft Windows™ 10, and we will consider an enhancement to include this scripting step in the FracPaQ GUI. Alternative routes to preparing a node file in the required format include using digitizer software (e.g., Engauge Digitizer available on GitHub (https://github.com/markummitchell/engauge-digitizer)). Example node files for each dataset used in this paper are provided with FracPaQ. Note that the (x,y) coordinates are in the reference frame of the image, and their default unit is therefore the pixel. The range of pixels is determined by the x and y coordinate range in the input node file.
The alternative input file type is the image file from which FracPaQ will detect fracture traces automatically. The user can specify an 8-bit binary image file name (e.g., in. JPG/.JPEG or. TIF/ .TIFF format) containing fracture traces. The freely available ImageJ tool works well for producing binary images from photographs (e.g. Fujii et al., 2007). FracPaQ then uses the Hough transform method to find co-linear patterns of pixels in the image and produce fracture traces (Kemeny and Post, 2003). Note that the Hough transform method, when operating on an input image file, is restricted to finding fractures with a single linear trace, so multi-segment fracture traces will be not be created. The user can adjust the Hough transform processing parameters using four text boxes on the GUI, which become enabled when the file type is 'image file'. Two example image files for input are provided with FracPaQ -Macduff. tif and Orkney. tif.
Note that, in general, node files of fracture traces traced by the user will contain more detail and be more accurate than an image file processed by the Hough transform. We have included the image file/Hough transform input method to allow a 'quick look' at a fracture pattern based on a digital image. Validation of user input is currently restricted to basic format and range checking on the values entered into the text fields on the GUI. Input file formats are not currently checked in any detail and the user is best served by supplying a correctly formatted node or image file.

Internal data structure and calculated quantities
FracPaQ can quantify the lengths and orientations of the constituent fractures, and the intensity, density, connectivity and permeability of the fracture pattern. These attributes and their units are listed in Table 1.
If the input file is valid, FracPaQ builds a MATLAB™ struct array of fracture traces (1 per fracture in the input file) composed of one or more segments delimited by nodes. The processing required to generate the selected maps and graphs operates solely on this hierarchical data structure of spatial coordinates. Viewed in this paradigm, the quantification of a fracture pattern becomes a problem in 2D coordinate geometry. Given the node coordinates of each trace and each segment, it is computationally trivial to calculate the lengths, orientations, and connectivity of the fracture network. Spatial densities, including intensity & density, are estimated using the circular scan window method of Mauldon et al. (2001). The permeability of the fracture network is estimated using a constant aperture, parallel plate model and the crack tensor approach (Oda, 1983;Suzuki et al., 1998). The data structure in FracPaQ contains the attributes necessary to calculate the 2D crack tensor.

Generating maps and plots
After selection of an input file, the user clicks the Preview button to load a fracture trace map into the main plot window, and populate the data structure. Statistics calculated from the data input file are displayed in the text box (lower left of the main window, shown as a portion of Fig. 2), including the x-and y-coordinate ranges, the number of traces, segments and nodes. Many graphics packages use a coordinate reference frame with the origin in the top left corner. The main plot window in FracPaQ and the default setting for the maps in the figure windows is to use an origin in the standard cartesian location of the bottom left corner. Data files produced in graphics packages such as Adobe Illus-trator™ are displayed with a 180 reflection about the X-axis when Previewed in FracPaQ. To correct for this artefact, the user can click the Flip Y-axis button to reset the coordinate space origin to the top left corner. Repeated clicking of the button simply toggles the direction of the Y-axis (normal or reverse). The current setting of the Y-axis displayed in the main FracPaQ window is the one used to produce any requested output maps.
The next step is to select output options on the right hand side of the main window, and then click Run to produce the maps and plots. Each selected map or plot is displayed in a separate figure window, and written to a separate. TIF file in the current folder. While each figure window is visible, the user can exploit the standard MATLAB functionality to resize or reformat the figure as they wish, and can save the figure to a filename or folder different from the default. The default length unit used in FracPaQ is the pixel. To produce maps and plots with lengths and coordinates measured in metres, the user enters a number in the text box "Scaling (pixels/metre)". That value is used in all subsequent maps and plots.
In the tests that we have conducted to date, run time has been very satisfactory, with most operations completed in a few seconds on standard desktop computers purchased within the last three years. The exceptions to this performance are the maps of estimated intensity and density, i.e. measures of spatial density, where run times are typically of the order of minutes rather than seconds. We have implemented a MATLAB™ WaitBar to provide basic progress information for these tasks, and we will focus future efforts on improving the performance of these calculations.

Approaches & methods
FracPaQ currently assumes that the input fracture traces lie on a statistically flat 2D surface, so that the effects of topography on the appearance of fracture traces does not require correction. The quantification of lengths and orientations is then reduced to simple operations in coordinate geometry. The graphs of length statistics are separated into plots for fracture trace lengths and fracture segment lengths. In FracPaQ, a fracture trace is a continuous line composed of one or more straight fracture segments. Fractures at any scale are believed to form through the interaction and coalescence of smaller fractures, so viewing the trace and segment data separately may be beneficial. Scaling of fracture lengths has attracted much analysis and debate (see reviews in Gillespie et al., 1993;Bonnet et al., 2001). Many researchers have performed a least squares fit to the log-log plot of lengths, and make the tacit assumption of an underlying power law (or 'fractal') distribution. This process is flawed for a number of reasons, and we prefer the use of Maximum Likelihood Estimators (MLE) to allow the user to assess the options for their data (see Clauset et al., 2009 for a discussion). The parameters controlling the most likely estimates for each type of distribution are written to the MATLAB Command Window. Fig. 3. A flow diagram for the FracPaQ application. The graphical user interface provides buttons that map to the processes drawn in ellipses. Operation of the software is simple: after picking an input file with the Browse … button, the user can Preview the results in the main application window. The user then selects output options, including the maps and graphs and their various parameters. Clicking on Run produces these outputs, each one in a separate figure. A copy of each figure is also saved as. TIF file in the current folder. The user can then Exit, or Browse … to select another input file.
The orientation distribution in a fracture pattern is important for unravelling the tectonic history of the rocks and in controls rock mass behaviour with respect to attributes such as permeability and strength. FracPaQ provides several different plots to assess the 2D orientation distribution in the pattern. Fracture angles in FracPaQ are defined as the angle of the fracture segment measured clockwise from the Y-axis, for the default assumption that the Y-axis is true North. This choice can be corrected for other trace map orientations if the user supplies an angle (in degrees) in the text box "Rotation of Y-axis from North". Rose plots or Rose diagrams are a popular way of visualising orientations in 2D, but the most common linear form of the plot is inherently biased (Nemec, 1988). We plot the rose with the area of each sector proportional to the frequency of orientations to avoid this bias. FracPaQ only plots angles of fracture segments and not fracture traces, as drawing a straight line between the two end nodes of a trace may not accurately capture the trace angle.
The spatial density (sensu lato) of fractures is known to vary as a function of distance from larger structures and is a critical attribute for assessing the transport properties of a rock mass. Maps of spatial density can provide insight into the processes of shear fracture growth from the interaction and coalescence of constituent microcracks or small fracture (Moore and Lockner, 1995). FracPaQ provides two measures of spatial density calculated from the input 2D fracture data (Table 1). Fracture intensity, labelled P21 by Dershowitz and Herda (1992), has units of m À1 and is defined as the total length of fracture in a given area (hence units of m/m 2 ¼ m À1 ). Fracture density, labelled P20 by Dershowitz and Herda (1992), has units of m À2 and is defined as the number of fractures per unit area. We estimate these measures from the data using the circular scan window method of Mauldon et al. (2001), applied to the coordinate geometry of the fracture trace and segment network. Mauldon et al. (2001) estimated fracture intensity as n/4r, where n is the number of fractures intersecting the perimeter of a circle of radius r, and fracture density as m/2pr 2 where m is the number of fractures terminating within a circle of radius r. FracPaQ generates a 2D grid of evenly spaced circular scan windows to fit within the fracture trace map area, where the scan circle diameter is defined as 0.99 of the grid spacing in x and y to avoid overlapping scan circles. The code then calculates the intersections (n) and terminations (m) of the fracture segments within these circles, and calculates the estimated intensity and estimated density values for the centre of each circle. This grid of values is then contoured using the standard MATLAB triangulation function to produce the maps of estimated fracture intensity (P21) and estimated fracture density (P20). The number of circles can be adjusted by the user in the text box "Number of scan circles" (NB. this is the number of scan circles in each of the x-and y-directions, thus the total number of circles is this number squared). Selecting large numbers (e.g. >30) of scan circles can result in long run-times (i.e. several minutes even on a powerful computer). The optimum number of scan circles depends on the specific attributes of the fracture pattern (see Rohrbaugh et al., 2002 for a detailed analysis). In particular, it must be remembered that the choice of scan circle size is related to the abundance and spatial distribution of fractures, and therefore the 'best' circle size will be different in different parts of the fracture trace map (Rohrbaugh et al., 2002). As a result, the estimates of fracture intensity and fracture density currently produced by FracPaQ have spatially varying errors, where a fracture pattern in heterogeneous. A user could customise the FracPaQ code to automatically run a systematic series of analyses, using different combinations of scan circle radius, number of scan circles and/or scan circle position. FracPaQ exploits the coordinate geometry stored in the data structure to make these potentially tedious and timeconsuming tasks faster, easier and more reliable (less user error), and makes it easier to evaluate the influence of scan circle attributes on the resultant estimates of fracture intensity and fracture density.
Maps of fracture patterns can provide important constraints about the ability of a rock mass to conduct fluids in the subsurface. FracPaQ provides two plots directly relevant to studies of fluid flow in fractured rock, the connectivity triangle and the 2D-permeability tensor ellipse. Manzocchi (2002) introduced the ternary plot of fracture connectivity with the 3 vertices of a triangle denoting I, Y and X nodes in the fracture network. Nodes are classified as 'I' (for isolated ends of traces), 'Y' (for branch points, splays or abutments) or 'X' (for cross-cutting intersections). More connected networks will plot towards the lower Y-X tie of this diagram, whereas less connected networks will plot towards the I apex. FracPaQ loops through the whole data structure of fracture traces and segments and finds the mutual intersections (using line-SegmentIntersect.m). The relative proportions of I, Y and X nodes are then calculated with respect to the total number of intersections found, and the connectivity triangle is plotted. FracPaQ also plots two 'contour' lines of connectivity, for C L ¼ 2.0 and 3.57 where C L is the number of connections (intersections) per line (or trace), after Sanderson and Nixon (2015). C L is defined as 4(N y þ N x )/(N i þ N y ) where N j refers to the number (not the proportion) of nodes of type j.
FracPaQ also provides an estimate of permeability in 2D using the cubic law, a parallel plate assumption and the crack-tensor formulation of Suzuki et al. (1998; see also Brown and Bruhn, 1998). The crack tensor (Oda, 1983) incorporates information about fracture sizes, orientations and spatial densities in a single Table 1 Fracture and pattern attributes quantified in FracPaQ.

Attribute Units Comments
Trace length Pixels or metres Segment length Pixels or metres Segment orientation

Degrees
Default is to assume Y-axis is North. A correction (positive or negative) can be applied to geographically register the data in the Rose plot.
Intensity, P21 Pixel À1 or metre À1 Density, P20 Pixel À2 or metre À2 Connectivity n/a I-Y-X ternary (triangle) plot from Manzocchi, 2002. Permeability Pixel 2 or metre 2 Depicted as a 2D ellipse with semi-axes oriented parallel to k 1 and k 2 . Scaling of the ellipse is √k in each direction, i.e. permeability in the direction of flow, not pressure gradient. Default is to assume Y-axis is North. A correction (positive or negative) can be applied to geographically register the data in the Rose plot.
measure. The crack tensor can be evaluated as a tensor of 2nd, 4th or higher rank, depending on the available data. We use a 2nd rank tensor approximation, and calculate the anisotropy of permeability in 2D. The crack tensor we use is calculated as P ij ¼ (p/4) *r *R 2 *T 3 * N ij where r is the density of fractures (number per unit area), R 2 is the mean of the squared lengths of fractures, T 3 is the mean of the cubed apertures of the fractures and N ij is the orientation matrix (e.g. Woodcock, 1977). Then, following Suzuki et al. (1998), the permeability is k ij ¼ (l/12) * (P kk * d ij e P ij ) where l is a factor between 0 and 1, and d is the Kronecker delta. The units of permeability calculated in FracPaQ default to pixels 2 , but if a scaling factor is entered in the 'Scaling (pixels/metre)' text box then the units are metres 2 . k ij describes a 2nd order permeability tensor for fluid flow through the fracture network. We plot the ellipse for permeability in the direction of flow, taking the ellipse axes as √k 1 and √k 2 (see Long et al., 1982). For permeability in the direction of pressure gradient the code can easily be changed to plot axes with lengths 1/ √k 1 and 1/√k 2 . By calculating fracture permeability in this tensorial form, it is relatively easy to add the permeability arising from any host rock or matrix porosity in a dual medium approach (e.g. Dershowitz and Miller, 1995). Note that FracPaQ currently assumes a constant aperture applied to all fractures in the network. Dealing with aperture variation in the fracture network to include more realistic aperture distributions is a priority for future releases. The factor l is an empirical constant that can be used to model the degree of connectivity in the fracture network, e.g. a fully connected network has l ¼ 1. Note also that surface roughness of fractures is not currently considered in the estimate of 2D permeability.

Results & applications
In this section, we describe the various outputs from FracPaQ exemplified with four fracture datasets collected from outcrop and laboratory data. Two outcrops have relatively planar rock pavements, and the other is from a gently undulating area of rocky coastline. The fourth example is from a polished thin section, imaged in the scanning electron microscope. These datasets collectively range over 5 orders of length scale, and include a range of lithologies (granite, limestone, and two kinds of sandstone). They also cover a wider range of tectonic settings, including late tectonic granite, folded limestone, faulted sandstone and jointed sandstone from a fold-thrust belt. In each case, fracture trace picking was performed in a graphics package at a constant resolution (or 'zoom'), rather than zooming in and out to find every fracture.

Length statistics
4.1.1. Example dataset e fractured granite at Souter Head, Scotland, UK We mapped fractures in the Souter Head granite, exposed on the coast about 6 km south of Aberdeen (Scotland, UK), using an Unmanned Aerial Vehicle (UAV, or drone) fitted with a digital camera (Fig. 4). This pink medium-to coarse-grained muscovite biotite granite intrudes Dalradian (Neoproterozoic) metasediments deformed in the Grampian Orogeny. The granite is surrounded by an explosion breccia and is hydrothermally altered. The intrusion has been identified as a sub-volcanic expression of the 'Newer Granite' suite by Porteous (1973), which as a whole were dated at 420-400 Ma (Stephenson and Gould, 1995;Trewin et al., 1987). A fracture network comprising 3 type of brittle structures cut the granite: joints, faults and quartz-filled veins. At least some faults post-date the explosion breccia, while some joints are probably related to the earlier cooling history of the intrusion.
We used a standard DJI Phantom™ 2 UAV with the FC200 camera, recording images at 14 megapixels, to construct a high- resolution photogrammetric survey, from which an orthophotomontage was derived. The UAV was flown in a regular grid pattern at a constant altitude of 30 m above local ground level, which undulates between 0 and 10 m above sea level at this locality. The digital photographs were processed in Pix4D Mapper™ software to build a photomontage of the whole area, including lens correction. This photomontage was then ortho-rectified with 12 ground control points acquired using a hand-held GPS device. The final orthomosaic (Fig. 4a) has a resolution of about 5 cm.

Length plots in FracPaQ
Quantifying the length distribution(s) in a fracture network is a key task, and we provide several different outputs to allow the user to visualise the raw data and assess the underlying distribution. The user can select histograms and log-log plots of length data. Separate figure plots are produced for fracture trace lengths and fracture segment lengths (Fig. 5). A dropdown list box is provided to control the binning of data for each histogram.
When selecting the "MLE analysis" checkbox in the FracPaQ GUI, the user will generate 6 extra plots of their length data, together with the most likely estimates of power law, exponential and lognormal distributions (3 plots for trace length data, 3 plots for segment length data) (see Fig. 6). At Souter Head, the lognormal distribution for the fractures in the granite is seen to be the most likely (94.56%), based on a Kolmogorov-Smirnov test (Rizzo and Healy, 2015). FracPaQ provides multiple complementary views of length data from a fracture trace dataset, and allows the user to assess the underlying distribution as a function of one of the three distributions.  Porteous, 1973). a-b) Histogram and log-log plot of fracture trace lengths from the area. Note that the trace lengths span well over two orders of magnitude. The red lines on the graphs mark the shortest and longest trace lengths found. c-d) Histogram and log-log plot of fracture segment lengths from the area. Segment lengths necessarily span a shorter range than the traces, but the range is still over two orders of magnitude. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Example dataset e fractured limestone at Spireslack open cast coal pit, Scotland, UK
A similar UAV mapping approach was used to build a photomontage of a fractured bedding plane in the McDonald (or Hosie) Limestone in the Spireslack opencast coal pit. Carboniferous coalbearing fluvio-deltaic rocks in this area have been worked since the 19th century, and form part of the main coalfields of the Midland Valley (Leslie et al., 2016). The McDonald Limestone is 2e3 m thick at this locality and is exposed in an 800 m long pavement in the back wall of the main void at the Spireslack pit. The limestone formation dips at 30e40 to the SE and is cut by numerous joints, veins and faults. The eastern end of this pavement was mapped using a DJI Phantom 2 drone flying at 30 m above the ground surface. This survey produced a photomontage covering approximately 400 m Â 30 m, which was lens corrected and orthorectified using 12 GPS ground control points in Agisoft Photoscan™ software (Fig. 7).

Orientation plots in FracPaQ
Fracture segment angles are shown from two sub-areas of the Spireslack dataset (Fig. 8). The user can change the 'petal' size of the rose plot using a dropdown list box of choices. The distribution of fracture angles from the area shown in Fig. 7b is shown in a rose plot and a histogram ( Fig. 8b and d, respectively).
The rose plots and the histograms confirm the presence of a fracture population as a function of orientation abundance (005 ) in the area with faults (Figs. 7b and 8b) compared to the area with only joints (Figs. 7c and 8a). Judging by the angles shown in the rose plot in Fig. 8b, a preliminary interpretation would be that the strike-slip faults formed in an apparently conjugate set oriented 150 and 005 . FracPaQ is ideal for making rapid comparisons of specific attributes of different fracture patterns, or sub-areas of the same pattern. The code could be modified to systematically sample many different sub-areas of a larger area, and compare the attributes of orientations, lengths, etc.

Intensity and density maps
4.3.1. Example dataset e experimentally deformed arenite from Clashach, Scotland, UK We mapped a fracture network in a cm-scale thin section of Hopeman Sandstone, a Permian age quartz arenite (Edwards et al., 1993;Farrell et al., 2014). The sample was taken from the Clashach quarry on the coast of Moray (Scotland), and a cylindrical core plug 20 mm in diameter and 50 mm long was prepared and deformed in a triaxial apparatus at EOST Strasbourg. The confining pressure was 25 MPa, and the pore fluid pressure was maintained a 10 MPa. The axial stress was increased to 290 MPa until the sample failed in shear. A polished thin section was produced from the core plug and imaged in the scanning electron microscope at Glasgow. A photomontage of backscattered electron (BSE) images of the whole thin section shows the through-going shear fracture transecting the field of view from top right to bottom left (Fig. 9a). Quartz and feldspar grains (middle grey and light grey, respectively) preserve many microcracks in the vicinity of the shear fracture (Fig. 9a). Fig. 10 shows intensity and density maps of the experimentally deformed sample of Hopeman Sandstone. The fracture trace map shows the traced microcracks and the through-going shear fracture, and the locations of the scan circles used in the calculations ( Fig. 10a; 25 each in the X-and Y-directions). The estimated fracture intensity map (Fig. 10b) shows high intensity (large length of fractures per unit area) in zones extending broadly parallel to the Yaxis e i.e. in the direction of s 1 (most compressive stress) applied to the sample. This pattern of fracture abundance is consistent with the through-going shear fracture forming from the interaction and coalescence of tensile microcracks aligned sub-parallel to the sample long axis and s 1 . The density map (Fig. 10c) shows 6 or more discrete clusters of high density (i.e. a high number of fractures per square metre) with the greatest values concentrated in the hanging wall of the through-going shear fracture. Note also the apparently even spacing of the high-density clusters measured parallel to the sample long (Ye) axis. Fig. 6. Fracture length statistics from the Souter Head granite, near Aberdeen (Porteous, 1973). It may be tempting to fit, or even assume, a power law scaling relationship from the straight-line portion of the segment length plot (Fig. 5d). Maximum Likelihood Estimation (MLE) provides a rational, quantitative basis for exploring alternative underlying models of distributions (Clauset et al., 2009;Rizzo and Healy, 2015). a) MLE analysis of a power law distribution, showing a poor fit overall to the whole dataset. b) MLE analysis of an exponential distribution, showing a reasonable fit over most of the range, but a poor fit at higher lengths. c) MLE analysis of a lognormal distribution showing a very good fit to the data. The observations coincide almost exactly with the estimation, obscuring the model line.  Fig. 11 (see Watkins et al., 2015a andWatkins et al., 2015b). The selected outcrops comprise arkosic sandstones of the Applecross Formation cropping out in the Achnashellach culmination of the Moine Thrust Zone. These rocks, deposited in the Neoproterozoic, were thrusted, folded and fractured in the Scandian orogeny (approximately 420 Ma; Watkins et al., 2015a). Circles of 1 m radius were chalked onto the outcrops and photographs were taken with a digital SLR camera at maximum resolution (18 megapixels). The photographs were ortho-rectified using ArcMap™, then scaled and oriented using Move™. All fracture traces were drawn on to these images digitally using the line tool in Move™. Tab-delimited ASCII text files of (x,y) coordinates for the fracture nodes were then exported from Move™ for use as input data in FracPaQ.

I-Y-X connectivity and permeability tensor in FracPaQ
For the Torridon Group sandstone example (Fig. 12), two outcrop datasets were analysed in FracPaQ to determine fracture connectivity and to estimate permeability anisotropy. The connectivity plots (Fig. 12a and c) are similar, with site SD11B having a slightly greater proportion of 'X' nodes and a lower proportion of 'I' nodes than SD11A, suggesting greater connectivity overall. The estimated permeability shows the k 1 azimuth for SD11A (Fig. 12b) is parallel to the long, low-density fracture set (Fig. 11c). For the SD11B, dataset the k 1 azimuth (Fig. 12d) is oriented parallel to the shorter, high-density fracture set. If the connectivity at site SD11A (Fig. 12a) was slightly greater, the shorter fracture set may have been far more influential in controlling the permeability, meaning

Limitations
FracPaQ currently assumes that the input fracture traces lie on a statistically flat 2D surface, where topography does not affect the appearance of any fracture trace trajectories. This requirement may Fig. 9. Details of Hopeman Sandstone dataset from a sample taken from Clashach Quarry, near Hopeman (east of Inverness, Scotland) and deformed in the laboratory at Strasbourg. a) A montage of back-scattered electron (BSE) images from taken from a scanning electron microscope (SEM). The grey regions are quartz (mid-grey) and feldspar (pale grey) grains, the black areas are porosity and microcracks. Porosity before deformation was about 11%. b) Fracture trace map prepared from the montage shown in a). Over 1600 fractures were picked from this small area, mostly concentrated around the diagonal through-going shear failure surface. Fig. 10. Estimated fracture intensity (P21) and estimated fracture density (P20) maps for the Hopeman Sandstone thin section. a) Fracture trace map showing the locations of the scan circles used to estimate intensity (m À1 ) and density (m À2 ) using the method of Mauldon et al. (2001). The area shown is just under 1 cm 2 b) Estimated intensity (P21) of microcracks. Note that the locus of maximum estimated intensity is not parallel to the through-going shear fracture (marked with a red line), and the asymmetry of areas of high intensity with respect to the footwall and hanging wall of the shear fracture. Note that the MATLAB triangulation and contouring algorithm can produce artefacts e.g. the contour map of estimated fracture intensity shows non-zero values around (3.5, 10), and yet, fracture traces are absent in this region on the trace maps (Fig. 10a). c) Estimated density (P20) of microcracks. Note the asymmetry of areas of high density with respect to the footwall and hanging wall of the shear fracture (marked with a red line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) be considered a severe limitation for areas with significant relief. Work is in progress to quantify this effect, by picking traces of known orientation from areas with differing topographic relief. We believe the discrepancy between the actual and the calculated 2D orientation distribution is small for areas where the relief is less than 10e20% of the maximum length of the area, and where the relief is undulose rather than significantly convex or concave over the whole area.
The acquisition of fracture datasets with a range of more than 3 orders of magnitude in their lengths is not common and has been non-trivial (see Bonnet et al. (2001) for a fuller discussion). However, there is no limit (other than computer memory) to the area that can be traced and then used as input for the FracPaQ toolbox. FracPaQ, used in combination with digital field acquisition methods such as UAV (drone) mapping, can quantify fracture networks where the (x,y) coordinate range is [ 10 3 m and the resolution of individual fracture lengths is < 10 À2 m. One of the key reasons for collecting fracture data over a wide range of length scales is to document possible scaling behaviours. FracPaQ provides Maximum Likelihood Estimators for power law, exponential and lognormal length distributions to help researchers quantify the distributions underlying their data.

Scope for further improvements
In addition to removing any bugs, we plan to implement improvements and extensions in future releases of the code. These currently include: Fig. 11. Details of the Torridonian Sandstone dataset taken from Achnashellach in the NW Highlands of Scotland, UK. a) Outcrop photograph of a fractured bedding plane in the Torridonian Sandstone, with a chalked circle of 1 m in radius (compass clinometer left of centre). This circular window is sub-area SD11A shown in c). b) Outcrop photograph of a fractured bedding plane in the Torridonian Sandstone, with a chalked circle of 1 m in radius (compass clinometer top of centre). This is sub-area SD11B shown in d). c) Fracture trace map of area SD11A derived from the high-resolution version of the photograph shown in a). Over 440 fractures were picked. d) Fracture trace map of area SD11B derived from the high-resolution version of the photograph shown in b). Over 650 fractures were picked.
improving the estimates of fracture intensity and fracture density through direct calculations of trace segment length within the scan circle areas; and through new GUI elements and underlying code to allow the user to better constrain the choice of scan circle size and location; adding a plot to illustrate the block size distribution showing the range of block sizes in between the fracture traces; processing of multi-colour input files to separate different fracture trace populations, and their statistics, based on the input line colour; better estimates of 2D permeability based on observed aperture distributions, rather than assuming a constant aperture, or coupling aperture to length (Vermilye and Scholz, 1995;Gudmundsson et al., 2001); adding branch and node statistics of 2D trace maps using the topological analysis of Sanderson and Nixon (2015); incorporating crack tensor estimates using 2nd and 4th rank approximations in 2D, and eventually 3D (Oda, 1983).
More generally, we are also working to ensure compatibility with GNU Octave, the open source platform compatible with MATLAB™, and to improve the run-time of the computationally heavier processing tasks (e.g. calculating fracture intensity and density using the scan line method). One overarching aim for releasing the code in open repositories is to generate feedback and  (Watkins et al., 2015a). Both areas are from Anticline 1 in Watkins et al. (2015a) and are 1 m in radius. a-b) In this area (SD11A), the connectivity is dominated by I and Y nodes, and the 2D permeability anisotropy (k 1 :k 2 ) of 4:1 is oriented with k 1 trending 140 . This trend is sub-parallel to the longest fractures in the pattern, but not the most numerous fractures (see inset rose plot). c-d) In this area (SD11B), the connectivity is broadly similar in terms of I:Y:X ratio, but the permeability anisotropy (k 1 :k 2 ¼ 3:1) is aligned parallel to the shorter but more numerous fractures (see inset rose plot). Note that the fractures counts are different for the two areas, with n ¼ 445 and n ¼ 660 for SD11A and SD11B, respectively. Also note that permeability is plotted in the direction of flow and the principal axes are scaled as √k 1 and √k 2 as described in Long et al. (1982). suggestions from the wider community e both researchers and developers e to help us improve the toolbox over time.

Summary
This paper describes a new software-based toolbox for quantifying fracture patterns in 2D, FracPaQ. The tools provide an objective and repeatable set of methods for analysing fracture data from a wide range of spatial scales, rock types and tectonic settings. The underlying methods used to quantify the fracture attributes are all published, but we have brought them together into a single toolbox. We have shown how FracPaQ can quantify fracture patterns across different scales (5 orders of magnitude), from thin sections (mm to cm scale), outcrops (m to dm scale), to regional maps (km scale). The graphs and maps of fracture pattern properties allow the user to systematically quantify statistical and spatial variations, and make comparisons between different areas and across scales.
The examples presented in this paper demonstrate how a geologist can use FracPaQ to explore the scaling laws underlying their fracture lengths and compare the orientation distributions between neighbouring sub-areas in an objective and repeatable manner. Robust methods of quantifying pattern attributes such as intensity, density and connectivity are also provided, and will provide insight into spatial variations and their possible correlation with distance from major structures. A basic 2D estimate of permeability anisotropy is provided to enable the investigation of bulk permeability as a function of different fracture pattern attributes, such as length distribution, orientation distribution or connectivity.
The toolbox has been developed in MATLAB™ and the source code is publicly available on GitHub™ and the Mathworks™ Fil-eExchange. The application, source code and sample input files have been made available in open repositories so that other developers and researchers will optimise and extend the functionality, and that "given enough eyeballs, all bugs are shallow" (Raymond, 2001).