Quantitative characterization of a naturally fractured reservoir analog using a hybrid lidar-gigapixel imaging approach

Understanding the orientation distribution and spatial configuration of natural fractures is important because these structural discontinuities significantly influence the behavior of many oil and gas reservoirs. As such, they impact fluid flow (e.g., Wilson et al., 2011b) and geomechanical state of the reservoirs (e.g., Heffer, 2012; Couples, 2013). Therefore, it is a common practice to include the contribution of fractures into static and dynamic reservoir models and simulations (e.g., Wilson et al., 2011a, 2011b; Bisdom et al., 2014). Recent advances in the 3D seismic imaging and analysis now allow the construction of geometrically accurate models with depositional and structural architectures constrained at resolutions of tens to hundreds of meters (Caers et al., 2001). However, most structural heterogeneities that can significantly impact reservoir behavior manifest at scales below the conventional subsurface imaging thresholds (∼20 m; Mrics et al., 2005). Unobservable structural heterogeneities such as faults, fractures, and compression structures (e.g., stylolites and compaction bands) introduce significant anisotropy within the rock mass resulting in permeability corridors or baffles and/or barriers, with hydraulic conductivities typically several orders of magnitude higher or lower than the surrounding rock mass (Aydin, 2000).

The use of outcrop analogs to generate geological conditioning data is a common strategy employed within reservoir modeling studies (e.g., Enge et al., 2007; Pranter et al., 2008). Scalable observations can be made from outcrops that are below conventional seismic thresholds, while providing spatially (especially laterally) continuous data (Jones et al., 2011). However, conventional methods of outcrop fracture characterization suffer from deficiencies in relation to their ability to efficiently capture sufficient information about fracture geometry, intensity, and orientation that can accurately represent the overall characteristics of the rock mass. A typical analysis of fractures in outcrop consists of collecting detailed observations manually using window mapping or scanline surveys (Priest, 1993). Manual survey techniques, although precise, are labor and time intensive, with resultant data sets restricted by the lack of spatial continuity offered by the sampling domain, due to inaccessibility of the upper reaches of most outcrops. Terrestrial laser scanning (TLS; also known as terrestrial lidar) is a proven close-range remote-sensing technique for outcrop studies, which may serve to enrich observations of a geological, geochemical, or geotechnical in nature, through the addition of a geospatial component (Enge et al., 2007; Buckley et al., 2008; Burton et al., 2011; Olariu et al., 2011; Hartzell et al., 2014). The inherent millimeter to submillimeter resolution and accuracy of laser scanning means that a wealth of high-quality data can be collected and analyzed within a relatively short period of time (Buckley et al., 2008; Seers and Hodgetts, 2014).

In the present work, we use terrestrial lidar characterization of rock discontinuities that encompasses fracture characteristics such as orientation, size, density, and spacing. We review existing methods, both conventional and digital outcrop based, to derive important fracture parameters from rock exposures. Building upon these studies, we utilize a combination of terrestrial lidar and georeferenced gigapixel photomosaics to map exposures of the Mississippian Boone Formation, an analog for reservoir objectives in the Arkoma Basin, USA.

The outcrop that forms the focus of this study is located at the War Eagle quarry, 6 km northeast of Huntsville, Arkansas on Highway 412 (W93°41.227′, N36°7.168′). Approximately 15 m of strata within the upper part of the Carboniferous Boone Formation is exposed within the quarry. The exposed strata belong to the lower Mississippian sequence and correspond to carbonate-ramp deposits with varying conditions of energy and depth (Handford and Manger, 1990). The range of depositional environments viewed in the War Eagle quarry has previously been interpreted as being deep-shelf margin to open-marine shallow-shelf edge settings (Liner, 1979). Lisle (1983) recognized that uppermost parts of the outcrop belong to Short Creek Oolite member (Fig. 1). On the other hand, virtually no oolites are present at the base of the outcrop, suggesting a more characteristic “Boone” style of deposition. Stratigraphically, the outcrop represents a typical Upper Boone succession with penecontemporaneous chert layers and nodules at the base and the Short Creek Member with massively bedded and persistent oolitic strata within the uppermost parts of the quarry wall exposures. These units are overlain unconformably by Chesterian Hindsville Formation, which is observed at the top of the eastern quarry wall. This unit contains cross-laminated, fine-grained sandstones and angular chert clasts, possibly derived from the underlying Boone Formation.

Three lithofacies can be identified from the War Eagle quarry outcrops. The contacts between the lithofacies are observed to be sharp where the entire quarry is arranged in coarsening upward cycles (Fig. 1B). The lower unit (termed as “layer 1”) is ∼12 m thick and composed of wackestone and packstone with crinozoan and bryozoan fragments scattered in a matrix of microspar (Lisle, 1983). Anastomosing bands and nodules of cherts and stylolites are abundant. The middle unit (termed as “layer 2”) is ∼5 m thick and composed of wackestone, mudstone, and carbonate breccia. Chert is nearly absent throughout the rest of the section. The upper unit (termed as “layer 3”) is ∼7 m thick and only exposed in the east wall of the quarry. It is composed of oolitic grainstones with common occurrence of cross stratification (see Fig. 2 for outcrop photos). This interval represents the periodic establishment of higher-energy environment within a low-energy carbonate shelf setting (Zachry, 1979).

The outcrop exposed in War Eagle quarry belongs stratigraphically to the upper Boone Formation. The Boone Formation and its lateral Osagean equivalents dip west toward the Oklahoma-Kansas boundary, where they form the historically targeted reservoir intervals in the Mississippi Lime hydrocarbon play. Potential reservoir facies observed in the outcrop are the grainstones, grain-rich packstones, and some sandstones. The observed vertical thickness of the grainstones in the quarry walls is ∼4 m. Laterally, these facies are traceable as far as the quarry exposure permits. Other observed sites in the vicinity had similar stratigraphic arrangements with similar ratios, which may be of significant volume. Deposited on a carbonate platform as sheets, subparallel to the shoreline, it is expected that lateral continuity in both strike and dip directions would be in the order of kilometers, similar to the modern grain-rich shoals (Fig. 3, Sivils, 2004). It is suggested that the Boone Formation formed in a shallow-dipping, carbonate-dominated depositional system that can serve as an outcrop analog for age-equivalent, hydrocarbon-producing sequences in the subsurface (Shelby, 1986). Cyclicity-related reservoir-seal lithofacies types can be used as a basis for reservoir comparison and modeling in the subsurface. These may include age-equivalent reservoirs in the midcontinent oil province such as those in Mississippi Lime hydrocarbon play (Keith and Zuppann, 1993, Fig. 4). In addition to the depositional setting and age equivalency, the Boone Formation possibly could potentially provide some insight with respect to the evolution of the subsurface stress field. Variations on the fracture orientations and densities observed on the outcrop and how they manifest within sedimentary units exhibiting different mechanical properties (e.g., Young’s modulus) could provide insights into the structural character of equivalent subsurface reservoirs.

The terrestrial lidar platform used in the study was a RIEGL VZ-400 terrestrial laser scanner with an operable range of ∼2–400 m, a ranging accuracy of 5 mm, and precision of 3 mm (RIEGL Laser Measurements Systems, 2011). A Nikon D800, 36.8-megapixel digital camera with a calibrated 20 mm lens was also mounted on the scanner, enabling the simultaneous capture of auto-referenced imagery. In addition, the laser scanner data were supplemented using off-scanner high-resolution panoramic images captured using a robotic camera stage (Gigapan EPIC Pro). A Nikon 3100 digital single-lens reflex camera equipped with a Sigma 500 mm variable-focal-length telephoto lens was used to acquire gigapixel panoramic images (Fig. 5). The range to imaging targets varied between 15 and 35 m, which resulted in ground sample distance between 0.4 and 0.9 mm in the object space. In comparison, the ground sample distance for the images captured with the TLS-mounted single-frame Nikon D800 camera setup varied between 3.7 and 8.6 mm.

Digital outcrop data at the War Eagle quarry were collected over two field days during Spring season under similar lighting conditions. To minimize occlusions (scanning shadow) and maximize draping accuracy, point-cloud and image data were acquired from multiple positions (11 and eight positions, respectively). The distance to the outcrop for all scanning and imaging positions ranged between 15 and 35 m. The two data sets consisted of a total of 11 scan positions, which were merged into a common coordinate system with the aid of a network of retro-reflective targets placed on the quarry walls. Furthermore, lidar-derived point clouds were aligned to the magnetic north with the aid of the scanner’s internal compass and clinometers, thus enabling measurements from the outcrop model to be compared directly with manually collected orientation data (internal compass accuracy is typically 1°; RIEGL Laser Measurements Systems, 2011). Additional reflectors were used to mark the manual scanline transects enabling the virtual scanlines to be accurately matched (e.g., Seers and Hodgetts, 2014). The resultant point cloud contained over 3 × 10⁸ points with an average sampling interval of 1 mm. Furthermore, a total of 3600 high-resolution panoramic images were simultaneously collected using the robotic camera stage, required for the generation of gigapixel texture maps. The use of exposures in the War Eagle quarry is motivated by the favorable geometries that the quarry walls exhibit in terms of digital outcrop data collection and analysis. For example, varying wall orientations within the quarry reduce Terzaghi discontinuity orientation bias (Terzaghi, 1965), and the planar, near-vertical and horizontal exposures limit data shadows and expose an extensive network of subvertical fractures.

The construction of outcrop models follows a similar approach presented by Buckley et al. (2008), Fabuel-Perez et al. (2010), and Rarity et al. (2014), with modifications to accommodate the registration of photomosaic imagery, which is explained in the following section. The workflow generally comprises four phases: (1) the merge and alignment of individual scan data from the two surveys; (2) cleaning, consolidation, and optimization of scan data; (3) registering external panoramic imagery to the scan data; and (4) triangulating and texture mapping digital outcrop surface meshes with referenced high-resolution photo-imagery.

Merging of point clouds was achieved using the Multi-Station Adjustment module within RiSCAN PRO software package (RIEGL Laser Measurements Systems, 2011). First, individual scans from the first survey were coarsely aligned using manually identified common tie points. Finally, all 11 scans were co-registered into a common coordinate frame using a global iterative closest point algorithm. The root mean square error of the final alignment between all 11 scans was ∼3 mm. Subsequent to the merger of scan data, operations were conducted to clean the non-outcrop artifacts (e.g., vegetation and boulders on the ground), consolidate all data to a single working data set, and optimize the resultant data set to reduce computational expenditure. This was achieved by resampling the point cloud in an effort to reduce the areas of excessively high and low point densities. Five-millimeter point spacing for the point clouds proved to be sufficient to resolve the discontinuity features of interest, with the mapped photomosaic imagery providing a medium through which millimeter to submillimeter fine-scale discontinuities could be identified.

Since registration of the gigapixel imagery to scan data was a significant effort within this study, the details of this procedure are given in the subsequent section. Following the registration of panoramic images, the final step of construction of a photorealistic outcrop model was triangulation and texture mapping with the georeferenced gigapixel images. A Delaunay triangulation method was used to generate a triangulated irregular network (TIN) from the cleaned and edited point-cloud data. An optimization procedure was then applied, eliminating the appearance of sliver tetrahedral (e.g., Seers and Hodgetts, 2014). Finally, the registered photomosaics were draped onto the polygonal mesh using texture mapping, creating a photorealistic reconstructed surface displaying dense RGB texture patches derived from the Gigapan robotic stage.

Because panoramic images were taken independently of the laser scanner, they require registration to the lidar data. Registration of an image requires establishing common points between the image and point-cloud data. Using these common points, external camera position and orientation with respect to the lidar coordinate system can be calculated. A network of control points, which included retro-reflective targets placed on the outcrop segments within reach and selected natural targets such as sharp corners of rocks, formed the basis of the registration. RIEGL’s RiSCAN PRO software package was used for the panoramic image registration task.

Most common photogrammetric packages employ corner-point projection when transforming 2D image coordinates to 3D lidar coordinates under the assumption of Cartesian image coordinates. However, the imaging geometry of the Gigapan platform can best be approximated by a spherical projection plane, in which the focal point does not change, while the viewing angle changes as the robotic arm rotates around a fixed axis. Microsoft’s Image Composite Editor (ICE) panoramic imaging software was used to stitch individual images. A spherical projection plane was used to minimize distortions inherent to panoramic imaging. The resultant images were treated as single-frame images where a pinhole camera model can be used.

The photogrammetry package used in this study (RiSCAN PRO) employs pinhole camera models similar to the one used in the Open Source Computer Vision Library, OpenCV (www.opencv.org). The calibration parameters defining a pinhole camera model consist of intrinsic and internal parameters. Intrinsic parameters refer to the camera sensor properties that are supplied by the manufacturer. Internal parameters are used to describe an ideal pinhole camera and consist of focal length, center of projection, and optical lens distortion coefficients. Focal length parameter is recorded in the image file, and center of projection was assumed to be the physical center of each image. This is a good approximation because of the spherical imaging geometry. Optical lens distortion was not modeled within the photogrammetry package because it was assumed to be removed during the stitching process. Using the established tie points and the simplified pinhole camera model, the position and orientation of the panoramic images were calculated. Qualitative confirmation of calculated camera position was done via captured lidar scans of the camera setup.

It should be noted that the registration of panoramic photomosaic imagery is imperfect in nature without applying the correct mathematical model for multicomponent stitched imagery with shifted pose (e.g., Schneider and Maas, 2006). However, the close proximity to targets and applying distortion removing algorithms to the imagery prior to stitching gave satisfactory registration results, on the order of a few mm to a few cm in the object space (Fig. 6).

An advantage of this method is that the panoramic photomosaics offer significantly greater pixel density on the outcrop surface compared to scanner-mounted camera images captured from equivalent ranges. Some fractures manifest themselves on the outcrop surface without topographic expression. These discontinuities are typically truncated from the analysis of a lidar-derived point cloud or tetrahedral mesh data set (Seers and Hodgetts, 2014). The combination of high-resolution panoramic imagery with lidar-derived geological surface reconstructions allowed the recognition of fractures with no discernable topographic expression (Fig. 7).

The data used in the analysis of this study consist of fracture attributes collected with compass and smartphone and extracted from TLS-derived surface reconstructions. In plan view, the quarry walls form a near-continuous rectangular perimeter (Fig. 8), meaning that subvertical fracture sets that are oriented parallel to one exposure should be manifested upon the adjacent wall, mitigating the effects of orientation bias. The quarry walls can be examined from almost a 360° range of viewing angles, thus unlocking the potential to access a 3D fracture data set. Rock discontinuities were measured by a traditional manual scanline method (Priest, 1993) with a geologic compass (such as Brunton) and a smartphone (such as Apple iPhone 5s with SMT 3-axis gyroscope and Bosch Sensortech BMA220 3-axis accelerometer; application used: GeoID by Sang Ho Lee, Engineering Geology and GIS Lab.), as well as a terrestrial laser scanner (manual guided and semi-automatic). First, discontinuities were measured along several scanlines with varying lengths between 20 and 60 m along the four quarry walls with a geologic compass and a smart phone. The procedure for all the scanline surveys followed broadly that of Priest (1993). With the intention of reducing orientation bias, seven scanline surveys were conducted on exposures with varying orientations (see Fig. 8 for locations of scanlines). Fractures exposed in plan view on the quarry floor were also recorded (e.g., “Floor” sampling domain in Fig. 8B). A total of 156 orientation measurements were recorded from all manual scanline surveys, including the plan-view exposures. The locations of the manual measurements were registered on the lidar scans in order to compare measurements between the digital outcrop-based approach and the manual surveys. It should be noted that a direct comparison of all fracture parameters between manually and virtually obtained data sets was not feasible within the scope of this study. The slow data acquisition rate inherent to manual surveys and limited number of fractures that were accessible at the base of the quarry restricted the comparison to be made only for the orientation statistics.

In the present work, the fracture orientation analysis from terrestrial lidar point clouds involved two approaches: manually guided and semi-automated. Manual analysis consisted of direct extraction of fracture attributes through supervised digitization of fracture planes. This was achieved by manual mapping of a fracture plane from the point-cloud data by selecting points that represent a rock discontinuity surface and calculation of normal vector of the best-fit plane through the selected points (Fig. 9). Dip and dip direction of the planes represented by these normal vectors are then calculated using Eigen analysis. High-resolution laser scans, colorized with the gigapixel photomosaics, guided the operator in interpreting the fractures. The entire digitization was carried out by a single interpreter to avoid variability due to influence of different operators (Scheiber et al., 2015). The manual digitization was limited to the outcrop areas with good exposure and lidar coverage from multiple scan positions. A homogeneous sampling throughout the quarry walls was aimed to be achieved; however, certain surficial artifacts such as vegetation and staining prevented reliable interpretation of discontinuities. Areas with such characteristics were neglected.

The second method for digital discontinuity orientation analysis was semi-automated. It was based on extracting quasi-planar discontinuity surfaces presumed to be fracture planes, bedding planes, and faults (e.g., Sturzenegger and Stead, 2009). This was achieved using an orientation tensor analysis (e.g., Woodcock, 1977; Fernández, 2005; García-Sellés et al., 2011; Laux and Henk, 2015), whereby the normal vector of the least-square plane for a group of neighboring vertices is approximated by the eigenvector corresponding to the minimum eigenvalue of their (Eigen) decomposed covariance matrix (see Fernández, 2005; Seers and Hodgetts, 2016a, for more in-depth discussion). Thus, the orientations of the normal vectors to the planar scanned surfaces were calculated via an Eigen analysis of the covariance matrix of a local point neighborhood (Metzger et al., 2009). After the orientations of the planar surfaces have been calculated, artifact non-discontinuity planar surfaces (i.e., those associated with quarry walls and bedding planes) were removed (Fig. 10). Best-fit planes representing the remaining points are produced and analyzed similar to the manual guided technique explained above. An unsupervised classification technique, k-means clustering, was then used to identify the orientation distributions of the main fracture sets.

Although point clouds or surface reconstructions enable fracture set distributions to be rapidly determined, the reliance upon point-cloud or surface reconstruction as a medium for fracture identification has a disadvantage. Some fractures that manifest on the outcrop surface as a linear trace without any topographic expression will not be detected upon the mesh or point-cloud data set (Seers and Hodgetts, 2016a). Therefore, co-planarity of the outcrop surface may not be diagnostic of the fractures, resulting in an underestimated view of the analyzed fracture network. To account for potential bias arising from the omission of such fracture traces, we used gigapixel textured meshes generated from Gigapan-captured photomosaics. The high-resolution nature of gigapixel texture maps allowed interpretation of healed fractures that generally do not have enough surface expression to be resolved by most TLS surveys. The interpretation of fracture traces from the textured surface reconstructions was manual in nature, with features digitized using CAD-based tools (e.g., ArcGIS ArcScene and RiSCAN PRO). Similar to manual fracture orientation identification, a single operator interpreted fracture traces.

Seven sampling domains were selected for detailed window mapping, required for the derivation of fracture abundance. Window sampling regions were selected based upon the quality of the exposure, outcrop orientation (to minimize orientation bias), and the localized occurrence of the lithological units present at the War Eagle quarry field site. With the intention of reducing data loss due to not interpreting all of the fractures, a lower cut-off length of 10 cm was used based on preliminary analysis of fracture trace data. Attention was given to make sure the cut-off size was distinguishable using the TLS-mounted camera images. The operator interpreted fracture traces bigger than the cut-off length within the selected areas. Fracture trace-length distributions, spacing, and intensity were then calculated using the digitized traces in these areas. A summary of the parameters calculated in this study and the methods of calculation are presented in Table 1.

Robust analytical techniques were employed to delineate fracture sets derived from both manual (i.e., with compass and smartphone) and virtual (i.e., extracted from digital outcrop models) orientation data sets. Initial efforts focused on qualitative assessment of field observations of large (>1 m) and systematic fractures, using the stereographic projections. In addition, K-means clustering algorithm was used to define fracture set membership. Statistical parameters of spherical distributions of fractures derived manually and virtually are compared to validate the suitability of digital data sets. Direction of vector sum ⁠, degree of preferred orientation (R%), shape parameter (⁠⁠), and cone of confidence were calculated for the corresponding fracture arrays (Table 2).

Fracture size exerts significant control over the gross network connectivity (Priest, 1993) and is a critical parameter for the generation of stochastic discrete fracture networks. In an attempt to estimate fracture size, fracture traces were identified by supervised digitization in selected sampling domains on the digital outcrop model (by window mapping). Using the previously obtained orientation statistics, fracture length distributions were investigated on a set by set basis. Probability plots were constructed for individual sets and used to assess the goodness of fit of measured distributions to empirically observed trace-length probability densities (log normal: e.g., Mayer et al., 2014; exponential: e.g., Robertson, 1970; hyperbolic: e.g., Segall and Pollard, 1983). In order to determine the significance of a trace-length population conforming to a specific distribution, the Kolmogorov-Smirnov (KS) test was performed. If the significance value of the test result is less than 0.05, the null hypothesis that the data comes from a specific distribution is rejected (assuming the optimally fitted distribution has the significance value of 1).

Different fracture abundance parameters can be obtained using different sampling spaces (i.e., 1D and 2D), discontinuity dimensionality (i.e., count, trace length, and surface area), and survey methodologies (i.e., scanline, window map, and circular scanline: Dershowitz and Einstein, 1988; Mauldon, 1998; Rohrbaugh et al., 2002). Much of the terminology was laid out by Dershowitz and Einstein (1988) and Dershowitz and Herda (1992), who defined a number of measures of fracture abundance by the P_xy system, where x denotes the dimension of the sampling region and y the dimension of the feature being measured. In the present study, we documented fracture abundance per fracture set as P₂₀ (fracture count over sampling area) and P₂₁ (total fracture length over sampling area). These measurements were taken using trace maps digitized within the window mapping regions described above (see example in Fig. 11).

Fracture spacing was calculated per fracture set using virtual horizontal scanlines placed on the selected sampling areas because fractures oriented subnormal to the constructed scanlines have greater spacing than the true value. To remove orientation bias, a correction based on Terzaghi (1965) was applied. Table 3 summarizes the parameters used in window mapping and virtual scanline surveys, which yielded the fracture intensity and spacing calculations.

Fracture orientations were measured by using manual scanline surveys, supervised digitization on digital outcrop data (virtual [supervised]), and semi-automated attribute extraction techniques (virtual [semi-automatic]). These methods yielded a total of 156, 532, and 5609 individual orientation measurements, respectively. Stereographic plots of poles to the planes (Fig. 12) indicate the presence of three major sets within both manually and virtually derived fracture populations—namely, north-south, northeast-southwest, and southeast-northwest. Qualitative confirmation of these population sets is provided by plan-view imagery of top and bottom of quarry walls (e.g., Fig. 2B), whereas (K-means) cluster analysis provided quantitative conformation, which indicated that 78% and 73% of total fracture population of manual and virtual data sets, respectively, fall within 10° (on horizontal plane) to a cluster center:

• (1) Set 1: N-S–oriented subvertical discontinuities, conforming to 46% of manually, 42% of virtually (supervised), and 40% of virtually (semi-automatic)–derived total fracture population.

• (2) Set 2: NE-SW–oriented subvertical discontinuities, conforming to 22% of manually, 24% of virtually (supervised), and 22% of virtually (semi-automatic)–derived total fracture population.

• (3) Set 3: SE-NW–oriented subvertical discontinuities, conforming to 32% of manually, 34% of virtually (supervised), and 37% of virtually (semi-automatic)–derived total fracture population.

The results of the statistical analysis of manual and virtual fracture orientations are summarized in Table 4. Calculated direction of vector sum of fracture populations indicates strong agreement between the manually and virtually derived data sets with a maximum angular deviation of only 6° (set 2) in the horizontal and 3° (set 2) on the vertical.Both measures of variance (⁠ and R%) indicate that set 1 has the least dispersion of all clusters, while set 2 has the greatest dispersion. In addition, measures of variance between the manual and virtual data sets reveal that all equivalent clusters from virtual data sets have greater dispersion than their manually derived equivalents. The 95% confidence levels indicate that the sample sets are unlikely to be randomly distributed.

A summary of trace length and density statistics for each fracture set is presented in Table 5 and Table 6, respectively. A total of 1294 fracture traces have been mapped on the digital outcrop model, where 27% of total fracture population were set 1, 44% were set 2, and 29% were set 3. Short but abundant set 2 fractures (mean length: 0.38 m) resulted in higher fracture intensity (both P₂₀ and P₂₁) values.

The analysis of fracture length distributions across the quarry indicates that the trace lengths of all recognized fracture sets approximate to a log-normal distribution (Fig. 13). While sets 1 and 3 fractures had similar average lengths (0.48 m and 0.44 m), set 2 fractures were generally shorter (mean length 0.38 m). Despite their short average lengths, set 2 fractures resulted in the highest P₂₀ and P₂₁ densities, whereas fracture sets 1 and 3 have lower and fairly similar density measurements.

Measured spacing generally overestimates the “true” spacing due to oblique intersection of fracture planes and scanlines resulting in apparent spacing. Terzaghi correction (Terzaghi, 1965) was applied to the spacing data set to calculate “true” spacing. Bisdom et al. (2014) explored the error associated with the common application of the Terzaghi method, concluding that it accurately calculated true spacing for angles between fracture poles and scanline less than 70°. We adopted this in our analysis and excluded the data from scanlines that were more oblique than this value. Table 7 presents the summary of the orientations used in Terzaghi correction.

Fracture networks generally exhibit strong spatial heterogeneity, frequently linked to faults, folds, stress fields, or lithological trends (Gauthier et al., 2000). These geological drivers control the large-scale flow behavior (Hardebol et al., 2015). The War Eagle quarry is an appropriate case to investigate the possibility of a fracture heterogeneity driven by the observed lithological trends. Here we present fracture heterogeneity as a function of fracture spacing variations using a large number of closely spaced (0.25 m) virtual horizontal scanlines placed on the digital surfaces of selected detail areas (e.g., Fig. 11). Spacing was calculated per fracture set by averaging the measurements from multiple scanlines of the same stratigraphic level. Each bar represents a spacing value (i.e., length of scanline/count) averaged from multiple scanlines of the same stratigraphic level. We applied Terzaghi correction to account for the sampling bias.

Uncorrected fracture spacing averaged 0.68 m, 0.85 m, and 0.97 m, whereas corrected spacing averaged 0.49 m, 0.79 m, and 0.88 m for layers 1, 2, and 3, respectively. Although the spacing data seem noisy (values vary locally between two adjacent virtual scanlines), there are some trends that can be seen. For example, the spacing of fracture sets 1 and 3 are similar throughout the section, whereas the spacing of fracture set 2 seems to be decreasing toward the upper stratigraphic intervals (Fig. 14).

Terrestrial lidar, coupled with referenced gigapixel photomosaics, provides an effective medium for fracture identification and analysis from rock exposures. The results of this study confirm previous work that suggests that the application of high-resolution discontinuity analysis is able to produce data sets compared to traditional hand-measured surveys (e.g., Becker et al., 2014; Seers and Hodgetts, 2014; Laux and Henk, 2015). We also confirmed that the virtual data sets, especially unsupervised, produced much noisier measurements.

Manually and virtually collected fracture data sets documented the occurrence of three fracture orientation sets: N-S, NE-SW, and SE-NW. Similar orientation distribution of fractures can be seen in observations made by Hudson et al. (2001), Braden and Ausbrooks (2003), Hudson and Murray (2003), Hudson and Turner (2007), and Chandler and Ausbrooks (2010, 2015) across the wider depocenter of the Boone Formation in northwest Arkansas. This may indicate the observations made from the quarry may be scalable up to several kilometers. Notably, the faults displayed in Figure 4 indicate the presence of three major orientations (N-S, E-W, and NE-SW) which are very similar to the orientations of the fracture sets presented in this study. This may indicate that the local stress fields that developed the observed fracture network in the quarry may be parallel or subparallel to the far-field stresses.

The Gigapan approach presented here could also be advantageous as a standalone exercise using structure for motion (SfM) photogrammetry as a complimentary 3D data set to TLS. In recent years, significant progress has been made in photogrammetry applications in fracture network characterization. Paired with the availability of open source codes to automatically manage image overlap from each position and applying a pixel-based automatic characterization, our hybrid Gigapan approach offers researchers much greater ground pixel coverage in their workflows (Hardebol and Bertotti, 2013; Zeeb et al., 2013a, 2013b; Healy et al., 2017).

The success of the virtual fracture characterization method largely depends on the scan and overlaid image resolution, as well as the size, geometry, and quality of the exposure. The minimum number of observation points to calculate the surface orientation with considerable accuracy depends on the scan resolution. Therefore, a prior knowledge of apparent fracture size distributions in outcrop would be advantageous to optimize the lidar survey for the extraction of rock discontinuities. It should be noted that the minimum fracture length that can be detected using the Gigapan setup is smaller than what can be detected using a single-frame TLS-mounted camera as demonstrated in Figure 7. However, low-contrast and/or minimal aperture fractures may still be elusive.

Fracture orientation measurements obtained from virtual and manual surveys revealed notable differences of dispersion. Calculated measures of variance (⁠/R%) indicate a higher dispersion for all of the virtual clusters compared to manual data sets. This may be a result of the high level of noise obtained within the virtual data sets. Although lidar offers rapid data acquisition rate, facilitating the collection of statistically significant data sets, the inherent dispersion due to noise with the digital outcrop-based techniques may impact the modeled network connectivity, resulting in underestimation of modeled network permeability (e.g., Berkowitz et al., 2000; Seers and Hodgetts, 2014).

Although noisier than the traditional hand surveys, orientation statistics of virtually derived fractures were in excellent agreement in terms of fracture strikes with those obtained manually suggesting that they can be potentially used as conditioning data for discrete fracture network models. Therefore, the virtual discontinuity analysis workflow presented here can be included as an alternative to the fracture modeling workflows. It should also be noted that the exposed fractures have their own variability, and therefore the virtual method may be accurately sampling the natural variability that the manual methods might have failed to detect. If we assume this is the case in the quarry and that a more accurate description of the fracture network can be obtained with the semi-automatic method, then this implies that a fracture set with similar strike and variable dip angles (Fig. 12) has high probability to have high fracture connectivity that may control the bulk permeability of the medium.

The ability to identify fractures from a digital outcrop model often depends heavily upon the fidelity and resolution of its surface display of RGB color (Vasuki et al., 2014; Seers and Hodgetts, 2016a). In the present work, high-quality photomosaic images were draped directly onto triangulated meshes, providing improved texture maps in terms of pixel density when compared to maps generated from on-scanner camera images. This is done to limit the effects of truncation on digital outcrop models. However, the referencing of large quantities of digital images is a considerably large exercise. The extra amount of time required to capture and process the optical remote-sensing data decreases the advantage in terms of data acquisition rate for lidar platforms within our hybrid workflow. A more practical solution may be to use calibrated images of outcrops (at the expense of measurement accuracy) or to adopt a robotic camera mount–based automated registration workflow such as the one presented by Cline et al. (2011). It should be noted that in the present work, the advantage of gigapixel textures is not exploited fully, given that 36-megapixel scanner images generally provide sufficient pixel density at 20 m for fracture traces to be readily identified. However, this hybrid lidar-photomosaic approach holds great promise for cases where lidar data acquisition is conducted at a considerable range (i.e., many tens to hundreds of meters; e.g., Hodgetts, 2013 and references therein). Supplementing mid- to long-range lidar data sets with gigapixel overlays, captured as multicomponent photomosaics using a telephoto lens–equipped digital single-lens reflex camera, could widen their analytical scope, enabling features such as fracture traces and sedimentary structures commonly aliased from on-scanner imagery to be discerned.

Despite its small extent, the War Eagle quarry fracture network displays a certain degree of heterogeneity. Because of the absence of an obvious tectonic control that could be responsible for such heterogeneity, we assume that the observed heterogeneity is due to variations in the mechanical properties of the exposed rock types. Of the three layers exposed in the quarry, layer 1 averaged the minimum spacing measurements (i.e., maximum P₁₀ density). We interpret this as being the result of layered and nodular chert occurrences characteristic to this layer. Figure 14 suggests that the low spacing in layer 1 was mainly due to set 2 fractures being more abundant in the lower quarry walls because other fracture sets were fairly constant throughout the vertical succession. Set 2 fractures are the short, mainly stratabound fractures that were observed to be primarily hosted in the cherty intervals. Another line of evidence came from the fact that the observed thickness of the chert layers (between 0.3 m and 0.4 m) and average trace length of the set 2 fractures (0.38 m) were fairly close, suggesting concentration of this family of fractures in chert-bearing rocks. Similarly, the calculated P₂₀ and P₂₁ intensity values signify the importance of set 2 in terms of potential influence of fluid flow. Being stratabound, set 2 may—together with bedding—provide an efficient connected pathway for fluids. Therefore, chert-bearing rocks may be the target interval for extraction of fluids in an analogous subsurface system.

The hybrid lidar and gigapixel texture mapping approach used in this study has been demonstrated to be a powerful medium for the characterization of naturally developed fracture systems as outcrop analog. Combining terrestrial lidar scanning and panoramic imaging leveraged enhanced data acquisition rate, increasing the statistical size of fracture attributes and extended pixel density, thus enabling those fractures to be identified that would otherwise be aliased from on-scanner camera imagery. In parallel to the virtual characterization, field-based conventional surveys were conducted. Identified discontinuity sets using digital techniques were found to be in excellent agreement with those obtained manually, although lidar-generated orientation clusters were more dispersed than their traditional counterparts. Conventionally and virtually obtained data sets suggested three major discontinuity orientations that are scalable across a wider depocenter than the studied outcrops. The hybrid workflow allowed the systematic investigation of fracture abundance in different depositional facies observed in the quarry. Considering the absence of a tectonic driving force, the fracture network heterogeneity observed in the War Eagle quarry suggests variable mechanical properties of the depositional units. The success of the fracture characterization using digital outcrop models depends largely on the scan and overlaid image resolution. While modern on-scanner digital cameras may provide sufficient pixel density for fracture identification at close ranges, our hybrid lidar-photomosaic approach holds great promise for cases where lidar data acquisition requires a considerable distance to targets.

We would like to thank Repsol USA for funding and permission to publish the results of this study. We are grateful to the members of the Remote Sensing Lab of University of Houston: Ünal Okyay, Virginia Alonso de Linaje, Diana Krupnik, Casey Snyder, Alvaro Berrocoso, Lei Sun, Daisy Huang, Darren Hauser, Preston Hartzell, and Craig Glennie for their help during the field work. We thank Lionel White of Geological and Historical Virtual Models, LLC for providing license for their ArcGIS extension, Geoanalysis Tools. We are grateful to Dr. Francesco Mazzarini and Dr. Enrique Gomez-Rivas for their valuable contribution during peer-review process.