Enhancing fracture-network characterization and discrete-fracture-network simulation with high-resolution surveys using unmanned aerial vehicles

A workflow is presented that integrates unmanned aerial vehicle (UAV) imagery with discrete fracture network (DFN) geometric characterization and quantification of fluid flow. The DFN analysis allows for reliable characterization and reproduction of the most relevant features of fracture networks, including: identification of orientation sets and their characteristics (mean orientation, dispersion, and prior probability); scale invariance in distributions of fracture length and spatial location/clustering; and the distribution of aperture values used to compute network-scale equivalent permeability. A two-dimensional DFN-generation approach honors field data by explicitly reproducing observed multi-scale fracture clustering using a multiplicative cascade process and power law distribution of fracture length. The influence of aperture on network-scale equivalent permeability is investigated using comparisons between a sublinear aperture-to-length relationship and constant aperture. To assess the applicability of the developed methodology, DFN flow simulations are calibrated to pumping test data. Results suggest that even at small scales, UAV surveys capture the essential geometrical properties required for fluid flow characterization. Both the constant and sublinear aperture scaling approaches provide good matches to the pumping test results with only minimal calibration, indicating that the reproduced networks sufficiently capture the geometric and connectivity properties characteristic of the granitic rocks at the study site. The sublinear aperture scaling case honors the directions of dominant fractures that play a critical role in connecting fracture clusters and provides a realistic representation of network permeability.


Introduction
Fluid flow and contaminant transport through fractured rocks are increasingly studied using discrete fracture network (DFN) models. The strength of the DFN approach lies in the explicit representation and inclusion of the geometrical properties of individual fractures which allows for a better understanding of fractures contribution to flow and transport (Smith and Schwartz 1984;Renshaw 1999;de Dreuzy et al. 2001;Neuman 2005;Reeves et al. 2008Reeves et al. , 2013. Monte Carlo analyses are commonly used to address the spatial variability in fracture network properties and subsequent uncertainty in flow and transport behavior (Berkowitz 2002;Neuman et al. 2007;Reeves et al. 2008;Liu et al. 2016a). Different techniques can be used to generate DFNs, although most studies rely on fracture data collected from outcrops and/or boreholes and then utilize stochastic methods to reproduce site-specific fracture attributes (Andersson and Dverstorp 1987;Patriarche et al. 2007;Reeves et al. 2014;Follin and Hartley 2014). The outcrop-based model relies on a true representation of the fracture network as observed on the outcrop, and preserves the topology and geometry of the fracture network (Odling 1997). However, limitations of outcrop-exposurebased models and the complex three-dimensional (3D) nature of fractures make it impossible to fully characterize fracture networks (Neuman 2005;Reeves et al. 2012).
Stochastic approaches assume that any observed fracture network can be considered as a single statistically equivalent realization of a random and stationary process of network generation (Cacas et al. 1990). Fractures are idealized as straight lines or planar discs/polygons (Baecher et al. 1977), and their geometrical properties (length, orientation, aperture) are considered independent random variables following certain probability distributions (Munier 2004;Lei et al. 2017). Based on these assumptions, synthetic networks can be generated from statistical parameters derived from field measurements. Large data sets are needed to obtain reliable statistical parameters; for instance, fracture length commonly encompasses at least 2-3 orders of magnitude at field sites (Bonnet et al. 2001;Clauset et al. 2009). In practice, it is often challenging to obtain such a dataset due to scale limitations imposed by exposure constraints and the resolution of the survey equipment and/or methodology. As noted by Odling (1997), any fracture map represents a subset of the fracture trace population with a lower limit defined by the resolution of the mapping technique, and the upper limit defined by the size of the observation domain. These limitations naturally lead to data censoring.
Recent advances in the use of unmanned aerial vehicles (UAVs; commonly called 'drones') and image processing (structure from motion) offer the possibility of highresolution mapping of fractures, and thus detailed characterization of fracture networks (Vasuki et al. 2014). Traditional techniques of field mapping such as airborne or satellite imaging techniques, are often resource intensive, difficult to employ, and suffer from coarse resolution that limits detailed and inclusive characterization of smaller background fractures. Conversely, millimeter-to centimeter-scale resolution can be achieved by flying UAVs at lower altitude. Furthermore, kilometer scale surveys can be performed with UAVs, and can easily include areas where accessibility is a challenge. These capabilities of UAVs greatly improve the ability of researchers to collect essential fracture network data (Carrivick et al. 2013;Vollgger and Cruden 2016).
Natural networks often exhibit complex patterns of spatial clustering (Gillespie et al. 1993; Sanderson and Peacock 2019); however, most DFN studies do not take spatial clustering into account. This limitation arises from conventional outcrop maps that do not offer the resolution necessary to detect and quantify clustering patterns. In the absence of high-quality two-dimensional (2D) data to determine if fracture clustering is present, fractures are typically assumed unclustered and a uniform distribution of fracture centers is assumed (Baecher et al. 1977;Dershowitz et al. 1991). This implies that fracture location can be modeled as a Poisson process (Ross 2010). The Levy-Lee method or the multiplicative cascade process (MCP) can be used to replicate fracture clustering. In the Levy-Lee method, distances between fracture centers are randomly drawn from a power law distribution, with the exponent of the power law representing the fractal mass dimension (Clemo and Smith 1997). The Levy-Lee method is sensitive to domain size and the input and output fractal mass dimensions are often inconsistent for small domain sizes. The second approach uses a multiplicative cascade process to generate a probability field map from which fracture centers are drawn (Schertzer and Lovejoy 1987;Meakin 1991). Unlike the Levy-Lee approach, the MCP approach uses a correlation dimension which better captures geometric patterns of spatial clustering (Cowie et al. 1993).
Fracture aperture is an important geometrical property governing fluid flow at the single fracture scale but is the least studied at the network-scale. This stems from difficulties in obtaining reliable aperture measurements under in-situ stress conditions, as measurements at the outcrop surface cannot be directly correlated to subsurface apertures due to weathering and lack of confining stress, both of which lead to overly dilated fractures. Additionally, roughness along fracture walls makes it challenging to determine if mechanical apertures measured in the field are representative of fluid flow processes. Mechanical aperture is defined as the physical distance normal to fracture walls, whereas hydraulic aperture refers to the cubic law aperture for a fracture that is computed from a measured flow rate (Snow 1965). For most cases, roughness along the fracture walls leads to mechanical apertures that are considerably larger than hydraulic apertures (Raven and Gale 1985;Renshaw 1995;Chen et al. 2000). Several mathematical and empirical expressions have been proposed to correlate mechanical and hydraulic apertures (e.g. Witherspoon et al. 1980;Hakami 1995;Olsson and Barton 2001).
Several studies have focused on fracture characterization using UAVs (Hardebol and Bertotti 2013;Tavani et al. 2014;Duelis Viana et al. 2016), spatial clustering representation (Darcel et al. 2003;de Dreuzy et al. 2004), and fractures aperture modeling (Olsson and Barton 2001;Olson 2003); but combining all of these components remains challenging. For instance, Bisdom et al. (2017) combined detailed UAV mapping with sublinear aperture to length scaling; however, they used a deterministic approach to generate the DFN and the equivalent network permeability was not compared with field data.
This study investigates the use of UAVs to improve fracture network characterization and integrates these results into stochastic DFN simulations. The proposed integrated approach uses a multiplicative cascade process to reproduce fracture clustering, and models aperture according to a sublinear aperture to length relation. To assess the relevance of this approach, network-scale equivalent permeabilities obtained from stochastically generated DFNs are compared with pumping test results. Results of this study will serve as a basis to develop a watershed-scale surface water/groundwater model that investigates the impacts of climate change on water resources in northern Togo. Like other regions in West Africa, the subsurface in northern Togo is underlain by fractured rock aquifers with low storage capacity. Changes in land use, demographic stress, and projected climate variability over West Africa raise questions regarding groundwater sustainability. Specifically, the DFN results will be used in the spatial discretization of the model and to assign permeability values that capture the anisotropy and permeability of the fractured bedrock.

Study site
The survey site is in northern Togo situated between longitudes 0°07′-0°19′ and latitudes 10°48′-10°58′ (Fig. 1). The geology of northern Togo is made up of two domains: the West African Craton (WAC) to the north, and the Volta Basin to the south (Eglinger et al. 2017). The WAC consists of Archean to Paleoproterozoic crystalline rocks that were last deformed~2,000 Ma ago during the Eburnean orogeny (Kalsbeek et al. 2008). The Volta basin lies unconformably on the southeastern margin of the WAC and dips 1-10°to the southeast. Strata of the Volta Basin have been affected by the Pan-African orogeny (~600 Ma) and show an ENE bedding orientation. The Volta Basin strata are divided into three supergroups, which presented from youngest to oldest are: the Tamale, Oti, and Bombouaka supergroups. Only the Bombouaka supergroup has surface expressions in the study region. The Bombouaka supergroup consists of coarsegrained quartzitic and feldspathic sandstones, interbedded with clayey-silty layers. It is overlain by the Oti supergroup, which consists of shales and greywackes. The Tamale supergoup is largely composed of conglomerates.
Deformation features in northern Togo are the result of the Eburnean and Pan-African Orogenies. The Eburnean orogeny evolved in two stages (Feybesse et al. 1990). The first stage is a NW-SE crustal shortening (2.2~2.19 Ga) that generated Indicators of the transcurrent regime consist of NE-SW strike slip faults, NE-SW and NW-SE en echelon faults, and NW-SE tension veins (Lompo 2010). The Pan-African orogeny, unlike the Eburnean orogeny, was predominantly a compressional event with at least four phases (NE-SW, ENE-WSW, SE-NW, and SSE-NNW shortening); these are recorded by west verging thrust planes and conjugated strike slip faults (Tairou et al. 2012). The survey site is located in the WAC and consists of lightcolored granodiorite composed mainly of plagioclase, quartz, and amphibole. The texture is phaneritic and quartz veins are observed at the site (Fig. 2). The outcrop is highly fractured with no shear displacement observable. The joints (Fig. 2) have a primary NE-SW orientation and an irregular spacing pattern. The secondary fracture orientation is ESE-WNW. The crosscutting and abutting relationships observed ( Fig. 2) suggest that the ESE-WNW set formed prior to the NE-SW set. This location was chosen because it represents one of the places where the WAC crops out with a good exposure in northern Togo.
Given the low primary rock permeability, groundwater flow in the region is controlled by secondary permeability associated with fracturing and weathering. Results from single-well pumping tests yield a permeability range of 7.9 × 10 −15 -1.6 × 10 −12 m 2 in the crystalline rocks of the WAC (Appendix Table 5). The region has a Tropical Savanna climate (Köppen classification) with an average annual precipitation (1992-2016) of 1104 mm/year. The vegetation consists of sparse deciduous shrubs and woodlands which is ideal for outcrop surveys.

Methods
A drone survey was conducted to map fractures in well exposed, horizontal outcrop in the WAC. Fracture traces on the outcrop map were first digitized, and then followed by statistical analyses of fracture attributes (i.e., orientation, length, and spatial correlation). Fracture statistics were then used to generate multiple 2D synthetic network realizations. Fluid flow was computed within each network by solving head at fracture intersections. Two scenarios were used to model aperture: (1) constant aperture and (2) a sublinear aperture to length scaling. Principal components of the permeability tensor were averaged over all flow simulations to determine the network-scale equivalent permeability. The equivalent permeability was calibrated to single-well pumping test results by adjusting either aperture in the constant aperture case or the Young's modulus in the aperturelength scaling relation.

Drone survey
The study outcrop was surveyed using a DJI Phantom 4 Pro quadrocopter. The latter was equipped with a camera that features a 2.54 mm, 20 megapixels CMOS sensor. Table 1 gives the camera specifications. Prior to the survey, the Drone Deploy application was used to define flight grids. The side and front overlaps were set to 70% to ensure a seamless stitching of adjacent photos and an accurate photogrammetric model. Once the side and front overlaps are specified, Drone Deploy determines the flight speed and frequency at which photos are taken. Flight altitude is adjusted until the target ground resolution is achieved. The following equation gives the relationship between resolution and flight altitude (Vollgger and Cruden 2016): where GR is the ground resolution in m/pixel, H is the vertical distance above ground level in m, W s is the width of the camera sensor in mm, W i is the image width in pixels, and f is the focal length of the sensor in mm. In practice, before surveying, the outcrop is inspected to determine the smallest feature of interest. This information is then used to estimate the appropriate flight altitude. Flight altitude was set to 20 m which corresponds to a ground resolution of 0.63 cm/pixel. Camera locations during the survey were determined using the Phantom 4 Pro v2.0 onboard GPS receiver. The field site is located in a remote region of Togo without survey benchmarks to constrain the georeferenced grid with ground control points as commonly employed in UAV surveys (James and Robson 2012; Bemis et al. 2014;Vollgger and Cruden 2016). As an alternative, the orientation and length of seven targets were recorded during the survey and later used to assess the accuracy of the drone map. This approach allows one to derive error estimates specific to fracture length and orientation. In addition to the drone survey, manual outcrop measurements of fracture length, orientation, spacing, and aperture are recorded for verification. The measurement transects are set perpendicular to the main orientation sets for reliable estimation of true fracture spacing and equal representation of fracture sets. Length and orientation are measured for each fracture intersecting the transect.

Drone data processing
Drone images are processed using the photogrammetric approach (Vollgger and Cruden 2016) which uses overlapping photos to reconstruct 3D structures, orthomosaic, and digital elevation models (DEM). Agisoft Metashape software (version 1.5.2) is used to process images. The first processing step identifies common features (tie points) in adjacent photos, using Structure from Motion algorithms (Fig. 3). Once tie points are detected, photos are aligned, and the geometry of the scene is reconstructed by determining camera positions and orientation.  A sparse point cloud (3.5 × 10 5 points) is generated as a result of the first step. The quality of the photo alignment can be inferred from the reprojection error, with low reprojection error values indicating good photo alignment. The reprojection error in this survey is 0.337 pixel which is equivalent to a 57-μm offset between the original and aligned photos. For reference, the USGS (2017) recommends 0.3 pixel as reprojection error goal for general UAV imaging purposes. Length error relative to the seven targets (Fig. 4) measured in the field varies between 10 and 20 cm and orientation error varies between 1 and 3°( Table 2). The orientation error is lower than the tolerance limit (3-5°) generally accepted for orientation measurements (Vasuki et al. 2014;Novakova and Pavlis 2017). Fracture characterization for DFN generation involves best-fit parameterization of both fracture orientation and length to probability distributions, and hence, these small error values are inconsequential. Next, the sparse point cloud is densified by extracting the position of all points in the aligned photos. The dense point cloud (2.5 × 10 8 points) is then converted into a meshed surface that serves as the basis for deriving DEM and orthomosaic maps. Fracture patterns are analyzed using FracPaQ v2.4 (Healy et al. 2017). Given the small size of the survey area, fracture traces are digitized manually from the image for better results in FracPaQ. Automated and semi-automated fracture tracing algorithms are generally time efficient, but are prone to artifacts and thus quality control and manual adjustments are required (Gustafsson 1994;Bisdom et al. 2017). Sampling biases often associated with manual interpretations are minimized by digitizing at the same level of magnification (×1.5). When fractures are curved, they are interpreted as a straight line joining the two endpoints (Fig. 5). Fracture mapping is a classical example of natural data censoring. Fracture tracing is limited only to areas with good exposure, and no extrapolations were made for regions covered by soil. Fracture traces are then saved as scalable vector graphics file which serves as input for FracPaQ. The default length unit in FracPaQ is pixel, the conversion to SI unit is based on a scale factor (9.6 mm/ pixel), which is expressed as the ratio of the pixel dimension of the input image and the known field dimension.

Statistical analysis
The survey site has a low topographic gradient and the study outcrop lacks vertical exposure necessary to determine fracture attributes in the vertical dimension. Additionally, borehole geophysics were not performed on any of the wells in the area, including the wells that generated the pumping test data. Thus, the statistical analysis is restricted to 2D fracture data. The statistical analysis approach first defines fracture set orientations within the fracture population, analyzes for the distribution of fracture length, and examines trends in fracture spacing. Directional data are analyzed to determine the orientation of each fracture set. Fracture length data exhibit linear trends on a log-log plot and are fit to a power law distribution. Trends in the observed fracture clustering are captured by determining the two-point correlation dimension. Spatial occupancy (a surrogate for fracture density) is quantified using the fractal box-counting method.
Orientation sets are identified by plotting fracture strike on a rose diagraam. For each orientation set, the mean orientation (θ m ) is computed following the method proposed by Mardia (1972). To implement Mardia's method, orientation data are first organized into bins with the corresponding frequencies.
The mean orientation (θ m ) is given by: where where n is the number of orientation elements in a set, f i is the frequency, θ i is the bin value, and θ m is the mean orientation. The dispersion κ of orientations around the mean is computed using a maximum likelihood approach (Fisher 1995;Berens 2009). Previous studies have found that fracture length spans multiple orders of magnitude and is best described by a Fig. 4 Survey site showing the seven control target locations power-law (Davy 1993;Marrett 1996;Renshaw 1996Renshaw , 1999Odling 1997;Bonnet et al. 2001). Power-law distributions of fracture length infer that a characteristic length in the fracture growth process is undefined (Sornette and Davy 1991). This is advantageous for outcrop studies where a power-law exponent can be determined locally from smaller-scale data and extrapolated to describe fracture length at subregional scales. Lognormal distributions of fracture length have also been proposed, but these distributions can arise from censoring of the largest fracture length values due to limited outcrop exposure or lithological layering (Odling et al. 1999). In practice, the inverse cumulative distribution of fracture length is plotted, also known as a survival function (Ross 2010). If length exhibits power law scaling, the tail of the largest values will exhibit a linear trend on a dual logarithmic plot. The slope of the linear trend represents the power-law exponent, and the minimum fracture length value can be determined from the lower threshold of the power law trend.
Accurate characterization of the spatial organization of fractures within a rock mass is important for reproducing fracture connectivity and quantification of fluid flow. In this study, the fracture network exhibits spatial clustering, and the computed correlation dimension is used to describe the spatial organization of fractures. The correlation dimension (D c ) is determined using the two-point correlation method. First, fracture centers (barycenters) are determined and the distance between pairs of fracture centers is computed Bour and Davy 1999). The two-point correlation function, C 2 (r), is then determined as follows (Hentschel and Procaccia 1983): where N p (r) is the number of pairs of fracture centers separated by a straight line distance less than r, and N is the total number of points (barycenters). By plotting the correlation function C 2 (r) against the radius r, on a dual logarithmic scale, the two-point correlation dimension (D c ) is given by: A classical box-counting method is used in this study to compute the fractal dimension to provide a scale-invariant estimate of fracture density . It is based on fitting the domain of interest with a square grid with box size r. The network is scanned and the total number of boxes (N b ) that contain at least one fracture within the domain of interest is determined (Barton and Larsen 1985). This process is repeated for different values of r. The fractal dimension (D) is obtained by plotting the number of square (N b ) against the box size r, on a logarithmic scale. The fractal dimension corresponds to the slope of the line: where a is the intercept and D is fractal dimension. The boxcounting method is applied to a subset of the outcrop where soil is not present and fractures are well exposed.  Klimczak et al. 2010). Given the heterogeneous nature of fracture networks, it is often challenging to establish a representative elementary volume (REV). A multiplicative cascade approach (MCP) is used to reproduce the complex spatial distribution of fractures observed in the network by assigning fracture locations. The MCP is an iterative process of fragmentation of a shape into nsubdomains (Schertzer and Lovejoy 1987;Darcel et al. 2003). The MCP is advantageous because it preserves the observed fractal spatial clustering. A detailed description of the MCP method can be found in Moein et al. (2019); a summary of the method is presented here. The number of iterations is a function of domain size and minimum fracture length: where l m is the minimum fracture length, L is the domain size, and m is the number of iterations. For the first iteration, the domain is subdivided into four squares which are assigned probabilities such that: where P i is the probability, τ is the ratio of sub-domain side length to domain side length, and D c is the two-point correlation dimension. For subsequent iterations, each initial square is further subdivided into four, with new probabilities satisfying Eqs. (9) and (10). The final probability is the product of the parent and daughter probabilities (Fig. 6). A probability density map generated from a total of eight iterations is used to reproduce the observed spatial clustering. Fracture centers are extracted from the probability density map using a discrete inverse transform method. To implement the discrete inverse transform method, the probability matrix of the previous step is reorganized into a vector and the cumulative sum of the vector is computed. Next, a vector U containing n uniform numbers between 0 and 1 is generated; n is the number of fractures desired. Each element U i of U, is matched to the closest element of the cumulative probability vector, based on the following inequality: where C is the cumulative sum probability vector. Some elements of U are likely to give the same element C i which implies an unrealistic duplication of fracture centers if assigned to the center of the cell. To avoid duplicates, the location of fracture centers is drawn from a uniform distribution over the selected cell. Once fracture locations are determined, fracture orientation, length, and aperture are assigned. Fracture orientations are modelled using the Von Mises-Fisher distribution: where I 0 (κ) is the modified Bessel function of order zero, μ is the mean direction, and κ is an inverse measure of the dispersion around the mean. For each orientation set previously identified, random variates are generated using the enveloperejection method (Best and Fisher 1979;Berens 2009). The relative abundance of each orientation set is maintained by honoring posterior probabilities. Fracture length is modeled as a power law (Davy 1993;Bonnet et al. 2001): where F(l) is the cumulative distribution function, l min is the lower power law cutoff and α the power law exponent. A vector U containing n random numbers between 0 and 1 is first generated; n is the number of fracture centers. Length variates are then computed as: where L is the length variate vector. Once all parameters are specified, fractures are added to the network until the observed fractal dimension (D) and density (P 21 ) are satisfied. The latter is given by: where A is the area of the domain and l i is the length. The last step consists of assigning aperture to each fracture. Two methods are used to assign aperture: (1) constant aperture applied to all fractures, or (2) aperture is assigned using a sublinear aperture to length scaling model (Olson 2003;Schultz et al. 2008). The sublinear aperture to length is used because it is only applicable to open-mode I fractures (joints), which is the case at this study site. Furthermore, the sublinear relationship is consistent with the linear elastic fracture mechanics and supported by a compilation of joints and veins data sets. The aperture to length scaling model relates fracture length L to maximum displacement by (Schultz et al. 2008): where D max is the maximum displacement orthogonal to fracture walls at the center of the fracture and η is a constant of the form: where K Ic is the fracture toughness in MPa.m 1/2 , ν is the Poisson's ratio, and E is the Young's modulus in GPa. In the absence of mechanical properties specific to the study site, estimates from other studies were used. The constant η is estimated using a K Ic of 1.545 MPa·m 1/2 (Atkinson 1984) and a Poisson's ratio of 0.26 (Tugume et al. 2012). Young's moduli for granite/granodiorite range between 69 and 91 GPa (Lanaro et al. 2006). These values are specific to granitic rocks characteristic of the study region.
Representing joint aperture as an idealized ellipse, Olson (2003) determined that average opening-mode displacement D avg is related to D max by: Substituing Eq. (16) in Eq. (18) gives: Assuming that fracture walls are modeled as two parallel plates (Snow 1965;Witherspoon et al. 1980), D avg in Eq. (19) is equal to hydraulic aperture.

Flow simulations
The flow simulations utilized the DFN model developed by Parashar and Reeves (2012). This model assumes fully saturated fractures (necessary for implementation of the cubic law) and steady state conditions. Furthermore, it is assumed that the matrix is impermeable and flow occurs only through fractures. First, the hydraulic backbone for each fracture network realization is determined using geometric and hydraulic methods that remove dead-ends fractures and isolated clusters (Fig. 7). This is accomplished by sequentially numbering endpoints of all fracture segments (a fracture segment is defined as the portion of a fracture spanning between its intersection with other fractures of the network), and then scanning the DFN domain to filter out fracture segments that have at least one end-point unconnected to any other internal node. The volumetric flow within individual fractures is given by Darcy's law where Q is the volumetric flow rate per unit fracture width normal to flow direction (m 2 /s), T is the transmissivity of the fracture in m 2 /s, and ∇h is the hydraulic gradient. Next, hydraulic head is computed at each internal node; internal nodes are defined as intersection points between fractures. Any given internal node can be connected to a maximum of four nodes; however, connections between nodes is highly irregular due to spatial clustering and can lead to sparse, irregular coefficient matrices . Equation (20) and the conservation of mass principle are applied to all internal nodes of each network. This results in a set of linear algebraic equations where the hydraulic head at the internal nodes is the unknown. These equations are solved iteratively using the minimal residual method (MINRES)-a variant of the conjugate gradient method that can be applied to symmetric indefinite systems-until an internal convergence criterion (relative residual error = 10 -7 ) is reached. Once hydraulic head is determined at each node, Eq. Two scenarios are used to simulate flow. In the first scenario, representative of the conceptualization of a constant value of aperture, all fractures are assigned the same transmissivity according to the cubic law (Snow 1965): Fig. 6 First and second iterations of the multiplicative cascade process. At each iteration, a parent cell is fragmented into four daughter cells. A total of eight iterations is performed for each network realization where T is the transmissivity in m 2 /s, ρ is water density in kg/ m 3 , g is the gravitational acceleration in m/s 2 , μ is water viscosity in kg.m −1 .s −1 , and b is the hydraulic aperture in m. In the second scenario, the sublinear aperture to length relationship is used to relate each fracture to its hydraulic aperture to compute a distribution of transmissivity values. For each network, two constant head boundaries and two linearly varying head boundaries are used to simulate flow through the networks (Fig. 8). These boundary conditions are used to ensure a constant hydraulic gradient distribution throughout the DFN domain. The directional components of the transmissivity tensor are computed by solving for Darcy's law first in the x-direction and then rotating the network 90°to solve it in the y-direction (Zhang et al. 1996): Flow values across each domain boundary are recorded for each boundary configuration per network realization. Given that DFN inputs (length and orientation) are drawn from probability distributions to account for variability within geometric and hydraulic properties, several simulations are needed to account for the variability in the simulated network-scale transmissivity. A total of 50 flow simulations are averaged to determine the equivalent representative permeability tensor.
The principal permeability components k xx' and k yy' are obtained by taking the eigenvalues of the 2D permeability tensor: where k xx' and k yy' are equivalent to the maximum (k max ) and the minimum (k min ) permeability, respectively. The principal direction (ω) of the permeability tensor is obtained as follows (Zhang et al. 1996): The average equivalent permeability is given by: The geometric mean of permeability values from 24 single well pumping tests serves as a calibration target for networkscale equivalent permeability values. Due to the high transmissivity and low storativity of fractured rock, the effective radius of drawdown, although unknown, most certainly provides permeability estimates from a significantly greater scale than simulated by the DFNs. Theoretical studies have shown that 2D simulations often underestimate permeability, compared to 3D simulations (Min et al. 2004;Lang et al. 2014). This discrepancy is the result of the inherent difference between 3D and 2D network connectivity. For field applications, however, fracture outcrops often lacks vertical exposure and thus, DFN characterization is reduced to 2D. In the absence of vertical exposure, constraining a 3D stochastic DFN becomes challenging. Permeability values from pumping tests were normalized by aquifer thickness to obtain a 2D equivalent permeability, prior to calibration. Furthermore, some studies in the WAC (Tshibubudze et al. 2009;Dabo and Aïfa 2010) indicate that fractures are mostly vertical to subvertical, which implies that connectivity in the horizontal plane controls fluid flow.
For the constant aperture case, a single aperture value is adjusted until the average equivalent permeability matches the pumping test results. The same approach is used for the sublinear aperture to length scaling, but in this case, the Young's modulus is the adjusted variable. The calibration is considered acceptable when the difference between the simulated values and calibration target is less than 1%.

Fracture orientation
A total of 1,184 fractures were traced and extracted from the outcrop map (Fig. 9). The rose diagram suggests two predominant orientation sets: N000-039°(NNE-SSW) and N080-130°, which account for 39% and 26% of the fracture population, respectively (Table 3). The remainder of the fracture population shows two sets: N040-079°(18%) and N131-179°(17%). The two predominant sets (N000-039°and N080-130°) are consistent with those found at other locations in the West African Craton (Egal et al. 2002;Hein et al. 2004;Feybesse et al. 2006;Tshibubudze et al. 2009;Dabo and Aïfa 2010). The similarity in these findings suggests that these fractures originated from the same regional tectonic stress event. The preponderance of NNE-SSW trending fractures (39%) indicates the influence of the Eburnean orogeny (Lompo 2010). The two main sets are near orthogonal and this greatly increases probability for fracture intersection in the horizontal plane and the likelihood of forming connected networks for conducting fluid. The formation of interconnected networks is also confirmed by pumping tests in the WAC where rocks of granitic composition have little to no primary porosity and permeability.

Fracture length and spatial patterns
Fracture length at the field site spans over two orders of magnitude, from 0.07 to 21.21 m. It is highly likely that fractures longer than 21 m exist but are censored by the limited exposure; the ratio of soil area relative to the total area of the rock mass is approximately 15%. Longer fractures are more likely to occur in orientation set 1 (N000-039°) as evidenced by the lengthweighted rose diagram (Fig. 9). A power law trend is only observable over slightly greater than one order of magnitude with a power law exponent α equal to 2.2 and lower cutoff l min of 1.2 m (Fig. 10). The two-point correlation dimension D c is 1.47 (Fig. 11) and falls at the lower end of values (1.3 ≤ D c ≤ 2) reported in the literature (Bonnet et al. 2001), suggesting a highly clustered network. The strong clustering combined with the predominance of small fractures indicates a strong dependence on longer fractures to connect otherwise isolated clusters of fractures with higher local densities. For instance, 37.5% of generated network realizations are below the percolation threshold, mainly due to the absence of long fractures to connect isolated clusters. This observation is consistent with the findings of de Dreuzy et al. (2004) that suggest when 2 ≤ α ≤ D c + 1, the connectivity of the network is sensitive to the occurrence of longer fractures to improve the chances of forming a percolating cluster. The fractal dimension D is 1.6 and the fracture density is 1.7 m/m 2 in the subregion of the domain measured (Fig. 11). After testing, it is found that a fractal dimension of 1.6 is equivalent to 1,600 fractures in a domain of 40 m × 40 m with a density of 2.2 m/m 2 . Note the differences in fracture density for given scales that necessitate the use of a scale-invariant, space-filling fractal dimension for upscaling.

Flow simulations
Average permeability computed from the pumping tests (2.2 × 10 −13 m 2 ) correlates to a constant aperture of 57 μm and a Young's modulus of 87 GPa for the sublinear aperture to length relationship (Table 4). The calibrated aperture of 57 μm is in general agreement with the value of 65 μm used in other studies (Min et al. 2004;Baghbanan and Jing 2008;Liu et al. 2016b). The correlated aperture (sublinear scaling) ranges from 17 to 147 μm with an arithmetic mean of 29 μm. Fracture length follows a power law with an exponent of 2.2 and thus, the sublinear scaling aperture distribution is biased towards small values, which explains the discrepancy between the two values. The Young's modulus is within the range of values (69-91 GPa) for granitic rocks (Lanaro et al. 2006). The maximum permeability k max is oriented N049°for the sublinear case and N055°for the constant aperture case. The sublinear scaling case aligns best with the predominant orientation observed on the rose diagram. Orientation set 1 accounts for 39% of all fractures and thus, is more likely to contain longer fractures given the power law distribution of fracture length. Correlating aperture with length leads to the propensity for longer fractures along fracture set 1 and concentrated fluid flow in that direction. This effect is less pronounced in the constant aperture case because all fractures are given equal weight to conduct flow, and the overall flow distribution within the network reflect the propagation of the boundary conditions into the network structure (Reeves et al. 2013;Trullenque et al. 2017). In both cases, the maximum permeability does not align perfectly with orientation set 1, because fractures in orientation sets 2 and 4 also contribute to flow.
The horizontal anisotropy, denoted as the k max /k min ratio in Table 4, is 2.8 for the constant aperture case and 2.1 for the sublinear scaling case. If all fractures have equal aperture, the anisotropy of the permeability tensor is solely dependent on network geometry. Conversely, in the sublinear scaling case, the effect of the network geometry on the horizontal anisotropy is reduced since long fractures dominate fluid flow, and these could occur in any direction (Fig. 12). Fracture transmissivity follows a power law distribution with an exponent of 1.48 that scales as two-thirds of the length exponent (Fig. 13). This relationship arises from the cubic relationship between transmissivity and aperture, and the square root relationship between aperture and length: where b is aperture, L is length, and T is transmissivity. Combining Eqs. (26) and (27) gives: The distribution of transmissivity (F T ) as a function of the distribution of length (F L ) is obtained by the following transformation: where t is the transmissivity variable and ψ is a function such that: If length follows a power law of the form l -α , Eq. (28) yields: Manual scanlines vs. outcrop map Thirty manual measurements (orientation and length) from ten (10) scanlines are compared with the UAV map data. Although the scanlines were scattered throughout the survey area, the manual measurements are likely not fully representative of the entire fracture network. In addition, it was challenging to maintain a consistent way of interpreting fractures in the field vs. the UAV map. Figure 14 shows that the scanline method is biased towards long fractures (shown as a shift in the plot towards larger values) which is consistent with previous studies suggesting that small fractures tend to be underrepresented in scanline measurements (Mauldon et al. 2001;Watkins et al. 2015). Due to the small size of the manual measurements and the subjectivity in scanline locations, only the main orientation sets were recorded (Fig. 14). It took less than 20 min to collect the drone imagery over a survey area of 95 m × 95 m, an area that is not feasible for characterization using manual scanlines. In addition, the one-dimensional (1D) nature of the scanline method precludes spatial pattern analyses (i.e., fractal and spatial correlation  Case k eq (m 2 ) k max (m 2 ) k min (m 2 ) k max ÷k min Sublinear scaling 2.2 × 10 −13 1.9 × 10 −13 9.3 × 10 −14 2.1 Constant aperture 2.2 × 10 −13 2.0 × 10 −13 7.2 × 10 −13 2.8 dimensions) which according to the results exert considerable control on fracture connectivity and fluid flow. Overall, the scanline approach served to validate the UAV results; however, its relevance in fracture characterization was limited and required significantly more time in the field to conduct.

Conclusion
This paper describes a workflow that combines detailed UAV mapping with stochastic DFN simulations. Spatial clustering and its scaling in the characterized network were honored through the use of a multiplicative cascade process. Fracture aperture was modeled using the sublinear scaling relationship and a constant aperture case. The UAV survey was able to resolve features at the centimeter scale. This enabled a detailed statistical analysis of fracture geometry attributes (length, orientation) as well as spatial properties (two-point correlation and fractal dimensions) which are typically unknown for most fractured rock masses as the analysis requires highquality planar data that 2D transects cannot provide. Stochastic methods were then used to generate multiple realizations of synthetic networks and to solve for the permeability tensor. The equivalent permeability values required only minimal calibration to match pumping test results by changing the aperture value in the constant aperture case and the Young's modulus in the sublinear scaling case. These results suggest that geometrical and a b c d Fig. 12 Flow through individual fracture segments illustrate differences between the constant aperture cases (a, c) and the sublinear scaling case (b, d) for two network pairs. Color intensity is proportional to the amount of flow, where dark blue segments conduct the most flow. Note that flow is concentrated in long fractures in the sublinear scaling case spatial properties of fractures provide a reasonable framework to understand fluid flow in complex fracture networks. Furthermore, the results suggest that a sublinear scaling relationship is an acceptable aperture model in the absence of direct aperture measurements, and more importantly, likely captures the high degree of variability in transmissivity within field-scale DFNs, and better reproduces the dominant directions and patterns of enhanced fluid flow. Complex spatial distributions of fluid flow are particularly important for solute transport investigations.